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ABSTRACT 


GARRETT, LLOYD BERNARD. An Implicit Finite Difference Solution to the 
Viscous Radiating Shock Layer With Strong Blowing. (Under the direction 
of GEORGE LOUIS SMITH and JOHN NOBLE PERKINS. ) 

An implicit finite difference scheme is developed for the fully 
coupled solution of the viscous radiating stagnation line equations, 
including strong blowing. Solutions are presented for both air injec- 
tion and carbon phenolic ablation products injection into air at 
conditions near the peak radiative heating point in an earth entry 
trajectory from interplanetary return missions. A detailed radiative 
transport code that accounts for the important radiative exchange 
processes for gaseous mixtures in local thermodynamic and chemical 

equilibrium is utilized in the study. 

Starting with minimum number of assumptions for the initially 
unknown parameters and profile distributions, convergent solutions to 
the full stagnation line equations are rapidly Obtained by a method of 
successive approximations. Damping of selected profiles is required 
to aid convergence of the massive blowing cases. It is shown that 
certain finite difference approximations to the governing differential 
equations stabilize and improve the solutions. 

The present study results indicate lower wall radiative heat 
fluxes for carbon phenolic ablation than predicted by previous 


investigators. 
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INTRODUCTION 


A blunt spacecraft entering planetary atmospheres at earth hyper- 
bolic speeds encounters intense radiative heating rates, particularly 
in the frontal or stagnation region. Acceptable interior temperatures 
are maintained by mass transfer cooling through the use of heat shields 
constructed of polymeric ablator materials (Walberg and Sullivan). 

Near peak heating altitudes in many entry trajectories, such 
as entry into the earth's atmosphere from a direct return manned Mars 
mission, the mass of gas injected into the flow field is an appreciable 
fraction of the mass of the oncoming flow (Chin) and is sufficient 
literally to blow the viscous boundary layer off the surface of the 
spacecraft. 

This condition, which is generally referred to as strong or massive 
blowing, has a physically destabilizing effect on the flow field that 
seriously impairs numerical solutions to the governing equations (Libby 
and Sepri). Libby (1970 ) describes the physics of the problem to be 
that of an inner region near the wall which is dominated by pressure and 
inertia forces where viscous effects are small, and a thin viscous outer 
region which adjusts to the edge conditions. For small blowing rates 
the viscous boundary layer is near the wall and the presence of the 
solid wall has a stabilizing effect on the flow. However, for the 
massive blowing problem, the gaseous inner layer may not adequately 
stabilize the flow. 

The problem of massive blowing where the outer viscous flow adjusts 
to the edge conditions has been studied extensively for nonradiating 
flows (cf.,e.g. , Kassoy, Libby (1962, 1970), Libby and Sepri, and 
Kabuta and Fernandez). However, there has been limited attention 
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directed to the solution when radiation is coupled into the problem. Most 
previous approaches either required excessive computer time for the numeri- 
cal solution or were approximate analyses which lead to questionable 
inputs of the thermodynamic properties required for the radiation compu- 
tations. Since radiative heat fluxes (in addition to being strong func- 
tions of temperature and density) can also be sensitive to the chemical 
species within the shock layer, some of which are strong radiation emitters 
and absorbers (Hoshizaki and Lasher), care should be taken in defining 
these quantities across the entire shock layer. A complete discussion of 
previous approaches to the viscous radiating shock layer with mass addi- 
tion is developed in the chapter on Review of the Literature. 

The purpose of this investigation is to develop an approach for the 
numerical solution to the coupled viscous radiating flow field along the 
stagnation line of a blunt body under both weak and massive blowing con- 
ditions. It is required to solve the governing Wavier- Stokes equations 
without making unnecessary simplifying assumptions to the equations and 
in the numerical solution that could result in inferior inputs for the 
radiative flux computations. Since the flow equations are coupled, 
iteration or some multiple pass procedure is required. Thus, the problem 
becomes one of efficient iteration procedures. 

An implicit finite difference method, which has previously been 
shown to be computationally efficient for chemical nonequilibrium studies 
without mass addition (Blottner, 1969); is developed for the solution to 
the proposed problem. The nonlinear governing differential equations are 
written in finite difference form at all nodal points within the shock 
layer, with boundary conditions specified at the wall and immediately 
behind the shock. The formulation results in a tridiagonal matrix system 
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of algebraic equations which is efficient for machine computation. The 
governing equations are solved "one at a time" in succession across the 
entire shock layer. The overall numerical solution technique to the two- 
point boundary problem is by iteration of the flow variables at each nodal 
point in the shock layer. Since the governing equations are solved one 
at a time, rather than concurrently at each nodal point (as in initial 
value forward integration techniques), then the nodal spacing requirement 
for the overall solution is not limited by the stability requirement of 
the most unstable equation. 

A section is devoted to the stability of the finite difference solu- 
tions to the governing equations. Particular attention is directed to the 
stability of the species continuity equation with binary diffusion using 
central differencing and a two-point windward differencing scheme that 
provides automatic damping of the profiles. 

Radiation computations are carried out using an existing radiative 
transport computer program, RATRAP, developed by Wilson (1967)* The 
program uses the tangent slab approximation (one-dimensional) that accounts 
for absorption and emission within a layer of arbitrary optical thickness 
and is for equilibrium gaseous mixtures of hydrogen, carbon, nitrogen, and 
oxygen. Pressure, enthalpy, and elemental composition profiles are com- 
puted by the viscous flow field solution as inputs to the radiation program. 
Radiation fluxes are evaluated at various nodal points within the shock 
layer to provide coupling with the flow field. The transport properties 
for equilibrium air (i.e., viscosity and reactive Prandtl number, see 
Hansen) are used in the computations to account for the energy transport 
due to binary diffusion. 

Numerical, results are presented for both air-to-air injection and 
for the injection of the ablation products of a carbon phenolic ablator 
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heat shield into an air stream for a range of blowing rates of interest. 
Additional computations are presented for constant density flows and for 
viscous nonradiating flows with air-to-air injection. 

The air-to-air injection cases include calculations for a 3.05-meter 
spherical nose radius body entering the earth's atmosphere at 15*24 kilom- 
eters per second velocity at 6l kilometers altitude with a blowing mass 
rate to free-stream mass flow rate of one-tenth. The results are compared 
with an initial value forward integration method, integral approaches, and 
an "exact" numerical solution for identical free-stream and wall boundary 
conditions. An assessment is made of some approximations contained in 


these analyses. 
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REVIEW OF LITERATURE 

The scope of the literature survey is restricted primarily to the 
flow within the shock layer near the stagnation region of a blunt reentry 
spacecraft. The emphasis is on viscous radiating flows with mass addi- 
tion which is typical for hypervelocity enters into earth's atmosphere 
from manned interplanetary missions. 

Some of the earliest predictions for the radiative heat flux to a 
blunt body entering the earth's atmosphere at earth parabolic velocities 
or greater indicated that the heat flux at the stagnation point was pro- 
portional to the velocity raised to the tenth power (Meyerott). In more 
recent years the predictions of the levels of the radiative fluxes have 
almost progressively decreased, while the complexity of the computations 
have systematically increased. An excellent review of the advances in 
radiating shock layer analyses is given by Anderson in his 1969 paper. 
Figure 1, which is taken from his paper, is a diagram of the various 
assumptions used in the analysis of radiating shock layers. The organi- 
zation of this chapter will follow this diagram, with a section devoted 
to the radiation transport models and one devoted to coupled and uncoupled 
flow field analyses including the more recent approaches to the radiation- 
induced massive blowing problem. Much of the literature review is con- 
densed from Anderson's survey, and interested readers are referred to his 
paper and his extensive list of references for a more lucid description 
of contributions to radiation flow field analyses prior to 1969. 

Radiation Transport Models 

The salient features of radiation models as they are applied for 
tractable solutions to radiating flow field problems are discussed in 




6 


this section. The evolution of the radiation models which are incorporated 

into high-temperature shock layer analyses is traced. 

In general, the radiative energy exchange in a gaseous medium is 

governed by integral equations which involve temperature, mass density, 

and the number density of the individual chemical species integrated over 

both the radiation frequency spectrum and three-dimensional physical space. 

Hirschf elder, Curtis, and Bird ( p. 721) describe the functional dependence 

of radiation on these properties in the following manner. 

In the presence of radiation, a molecule has a certain 
probability of either absorbing or emitting radiation of 
a frequency characteristic of some transition from one 
quantum state to another. Or a molecule in an excited 
state has a certain probability of spontaneously emitting 
radiation of a certain frequency. Each substance therefore 
has an absorption spectrum which can be expressed in terms 
of the coefficient of absorption which is a function of the 
frequency. Since the absorption spectrum depends upon the 
distribution of the molecules in their various quantum 
states, then the absorption coefficient depends upon the 
temperature of the substance. The spectral lines for iso- 
lated molecules have a natural width (due to spontaneous 
emission of radiation), and as the molecules are brought 
together these lines become broader (due to pressure 
broadening) and become displaced (due to the distortion of 
the molecules themselves). Thus the coefficient of absorp- 
tion depends upon the density of the system . . . 

Since fluid elements (actually the individual atomic and molecular parti- 
cles ) both emit and absorb radiation, then the radiation exchange for both 
mechanisms must be considered. At a given point emission is a function 
primarily of the conditions at the point, whereas absorption is dependent 
upon not only conditions at the point but also is a function of the 
thermodynamic properties and the frequency of radiation emission of all 
the surrounding fluid elements. Radiation at a given frequency travels 
a "photon or radiation mean free path" before being absorbed, thus absorp- 
tion of the radiation is dependent upon the physical distance between the 
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emission source and the absorbing particle. Consequently, radiation 
exchange within a given volume is a function of three-dimensional physical , 
space. 

The problem of the general three-dimensional radiation exchange where 
there are many different chemical species or particles at various energy 
levels is indeed formidable, and simplifying assumptions are required in 
order to obtain tractable solutions to most radiation exchange problems. 

For stagnation streamline analyses, the enclosure volume is simplified by 
the "tangent slab approximation," that is, radiation heat fluxes are com- 
puted assuming that a one-dimensional planar slab of gas is present within 
which conditions remain constant except across the slab in the direction 
normal to the slab. This approach is almost universally applied in 
radiating shock layer analyses. At the stagnation line the justification 
for the assumption is that for bodies of large radii surrounded by a 
relatively thin shock layer, conditions vary slowly in the radial, direc- 
tion, whereas the major gradients are normal to the body. This effect 
has been investigated by Kennet and Strack, by Koh, and by Hoshizaki and 
Lasher. The investigations indicate that the error introduced by the 
approximation should be less than 5 percent. 

The tangent slab approximation permits one to evaluate the divergence 
of the radiative heat flux that appears in the energy equation by evalua- 
tion of the gradient of the heat flux in the normal direction only and 
uncouples the stagnation line solution from the rest of the shock layer 
as far as radiation is concerned. In the present analysis the tangent 
slab approximation is used. 

For radiation computations the gas is treated either as transparent 
or self -absorbing. A transparent gas is one which emits radiation but 
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does not absorb any incident radiation from the surrounding fluid elements. 
According to Vincenti and Kruger,, this approximation is valid only when 
the gas is optically thin at all wave lengths. An optically thin gas is 
one in which the characteristic mean free path of a photon is much larger 
than the thickness of the shock layer. Anderson points out that in 
practice the transparent gas assumption is reasonable only for reentry con- 
ditions where radiation effects first become noticeable, such as 10 km/sec 
entry velocities. For velocities around 15 km/sec, the transparent assump- 
tion can overpredict the radiative flux by factors of 2 or more. 

A self-absorbing gas both emits and absorbs radiation; that is, a 
fluid element locally emits radiative energy as well as absorbs energy 
from the surrounding fluid elements. The gas is treated either as gray 
or non-gray. A gray gas includes gray self -absorption which is a func- 
tion of temperature and pressure, but the absorption coefficient is 
assumed not to be a function of the wave length or radiation frequency. 

It was first indicated by Olstad and later demonstrated by Hoshizaki 
and Wilson ( 1967 ) that the non-gray model is by far the more realistic 
model and is important for most high-velocity entry missions. In this 
model the absorption coefficient is considered to be a function of wave 
length as well as temperature and pressure. 

An example, given by Anderson, that indicates the reason for includ- 
ing non-gray self-absorption is the phenomenon associated with high- 
temperature air. Air will absorb radiation in the vacuum ultraviolet 
(short wave length) region but is relatively transparent in the infrared 
(long wave length) region. For example, there are five orders of magni- 
tude variation with frequency in the continuum absorption coefficient 
of air at 14,000° K and one atmosphere pressure. Olstad' s 
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results for the inviscid flow of air and Hosizaki and Wilson's (1967 ) 
results for the viscous flow of air showed that substantial reductions 
in the radiative heat flux to the body were obtained in going from gray 
to the non-gray self -absorption models. 

The non-gray gas model must be applied to obtain realistic estimates 
for the radiative heat flux to the body and this requires a detailed 
integration of the absorption coefficient over the frequency spectrum. 

For air, significant differences still exist not only in the spectral 
absorption coefficients for certain chemical species but also in the 
radiation models that are coded for computer solutions (Suttles, 1971 )• 
However, computations dealing strictly with radiative transport in air 
for typical earth entry conditions are approaching a firm basis and 
simplified radiation models can be developed for absorption coefficients 
as functions of wave length, such as Callis' three-step model, to speed 
up the time-consuming radiative transport computations. When one also 
considers injection of ablation products into the air stream, radiation 
transport modeling becomes complicated indeed because of the added dimen- 
sion due to the presence of foreign radiating chemical species (Anderson). 
In the interest of quantitative results, it appears better at present to 
generate spectral absorption coefficients from the spectral details for 
the individual chemical species than to attempt absorption coefficient 
modeling on the basis of a specified overall chemical composition. 

Existing radiation transport computer programs are available which 
consider ablation products in addition to air chemical species, such as 
RATRAP, developed by Hoshizaki and Wilson (see Wilson, 1967 ); SPECS, 
developed by Thomas; and RADICAL, developed by Nicolet. These programs, 
which perform radiation computations based on the spectral details of the 
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individ ual! chemical species, are time consuming for stagnation line 
analyses. However, usage of these types of programs are required since 
realistic analyses have indicated that the addition of ablation products 
to the flow field can reduce the heat fluxes to the body by factors of 
two below pure air flow (see, for example, Hoshizaki and Lasher, 

Coleman et al. , and Chin). As was the case for air, differences exist 
in the computer codes due largely to uncertainties in the absorption 
coefficients and the spectral modeling of certain chemical species. 

RATRAP, which is used in this analysis, is somewhat time consuming, but 
it contains the appropriate detail to be compatible with the rigorous 
flow field competition expected of the present analysis. 

It is appropriate at this point to begin the discussion on flow field 
analyses with radiation. The radiation transport models discussion is 
terminated upon noting that differences in the heat flux predictions at 
the body can be a result of not only differences in the radiation trans- 
port codes but can also be a result of the inaccuracies in the solutions 
to the flow equations. In particular, the concentration or number densi- 
ties of strongly absorbing or emitting ablation chemical species can be 
strongly affected by computed temperatures and densities within the shock 
layer. 

Flow Field Analyses With Radiation 

The appearance of the divergence of the radiative flux as a term 
in the general energy equation governing the flow of a radiating gas 
couples the flow field and the radiative transport analysis. For rela- 
tively low entry velocities, the radiative heat fluxes are small and 
consequently do not exert any significant influence on the flow field 
thermodynamic or flow variables. However, for hyperbolic entry velocities. 
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typical of manned return missions from Mars, the radiative heat fluxes 
are large, and thus radiative cooling (energy losses) with the shock 
layer must he considered. The net effect of radiative cooling is to 
reduce the radiative heat flux to the body because of a reduction in 
temperature and a subsequent increase in shock layer density from adia- 
batic conditions (see, for example, Wick and Hoshizaki and Wilson, 19&5)* 
This coupling between radiative transport and gas dynamics for this 
problem is commonly referred to as the coupled radiating shock layer 
problem. 

Howe and Viegas were the first to investigate the flow of a viscous, 
radiating, self-absorbing gas in the stagnation region including the 
effects of mass addition. Since they assumed a gray radiation model, the 
radiation flux results are not quantitative. However, they showed that 
similarity considerations could be applied in the stagnation region of 
the viscous shock layer where radiation is present, thus reducing the 
governing partial differential equations to ordinary differential equa- 
tions. Howe and Viegas used a Levy-Lees type of boundary-layer trans- 
formation, which involves the integral of the viscosity-density product 
over physical distance, and similarity considerations to reduce the 
Navier-Stokes equations to compressible Falkner-Skan boundary-layer 
equations (see Schlichting, pp. 354) for axisymmetric stagnation region 
flow. They solved the momentum equation numerically by an initial value 
forward integration method (Adams -Moulton predictor corrector, see 
Hildebrand). The energy and species continuity equations were solved 
by numerical evaluation of the exponential integrals that appear in the 
exact analytical solutions to the differential equations. The solution 
to the coupled system of equations was iterated by converging on the 
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value of the shear stress at the wall. In their fully viscous analysis 
the mass injection rates were restricted to moderate amounts which did 
not upset the stability of the boundary layer near the wall. The diffi- 
culty with the stronger blowing rates can be traced to numerical instabil- 
ity problems associated with evaluating the exponential integrals in the 
exact analytical solutions to the governing equation. Wilson (197O), who 
applied a similar approach as that of Howe and Viegas (with the notable 
exception being his momentum equation solution), describes the problem 
for the fully viscous shock layer with massive blowing to be one of numeri- 
cal precision required to take the differences between the large numbers 
which appear in the exponential integrals. Further aspects of this 
problem will be discussed in association with Wilson's 197 ° paper and in 
the Results and Discussion chapter. 

In 1965 Hoshizaki and Wilson (1965) developed their integral method 
for the solution to the coupled viscous radiating shock layer about a 
blunt body. Fifth and sixth order polynomials were used to express the 
velocity and total enthalpy profiles, respectively, across the shock layer. 
In addition to presenting results for the stagnation region, they also 
presented results around the body using a forward integration technique 
with a li mi ted number of iterations on the shock shape. The solutions 
were restricted to no blowing and a transparent radiating gas. 

Hoshizaki and Wilson (1967) extended their integral approach to 
include injection of ablation products into the boundary layer for moderate 
blowing rates. And, they improved their radiation transport model to 
include specular (non-gray) self-absorption. In the solution, fifth and 
second order profiles were assumed for the velocity and elemental mass 
fraction of the ablation products, respectively. With the specified 
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velocity and species distribution profiles, they were able to solve the 
energy equation by using a method similar to that used by Stoiith and 
Clutter for boundary layer solutions. In the Smith and Clutter solution, 
the energy equation is solved across the shock layer by a method of super' 
position of two independent energy equation solutions. The nonsimilar 
terms in the energy equations (derivatives in the tangential direction) 
were included in finite difference form. Hoshizaki and Wilson (1967 ) 
assumed a binaiy diffusion model and assumed that the air and ablation 
products did not react chemically. Air transport properties were used in 
the analysis. They observed that, since the entire shock layer was 
treated as viscous, they did not have to match frequency dependent radia- 
tive fluxes at the viscous-inviscid interface. 

In 1968 Hoshizaki and Lasher extended the integral approach to the 
massive blowing problem (up to 10 percent of free-stream mass flow). The 
integral method was applied to obtain an approximate solution to the 
momentum equation for a fifth order polynomial representation of the 
velocity profile. The species continuity and energy equations were 
solved by means of similarity transformations and numerical integration 
of the resulting exponential integrals. Their detailed analysis of the 
spectral absorption coefficients for 20 air and ablation product chemical 
species showed that the carbon atoms, which diffuse far into the shock 
layer, act as strong radiation, absorbers. 

Chin, in 1968, developed a numerical method for solving for the 
radiation coupled inviscid stagnation flow with mass addition. In this 
analysis it was assumed that no mixing occurred between the ablation 
products in the inner inviscid region and the air products in the outer 
inviscid region which results in a distinct interface (in terms of 


chemical species) between the two regions. He integrated the conserva- 
tion equations for the air layer from the shock wave to the interface 
and from the body to the interface. He iterated on the wall heat flux, 
the blowing rate, and the velocity and enthalpy profiles until he con- 
verged on the enthalpy distribution and heat flux to the wall. It is 
interesting to note that the blowing rate is an implicit part of his 
solution and, for the conditions of interest in this analysis (U co ~15km sec, 
altitude ~ 65 km), he calculated a blowing rate of 7.6 percent of free- 
stream mass flow rate. His results were for a spherical body of radius 
256 cm, constructed of carbon phenolic ablative material. The solutions 
to the inviscid governing equations converged very rapidly; however, since 
the viscous region (which typically occupies about 10 percent of the shock 
layer for the Reynolds numbers of interest) was neglected, then there was 
no mechanism provided for the diffusion of the strongly absorbing carbon 

atoms and ions into the air layer. 

The first numerical solution to the exact Navier-Stokes equations 
for the thin shock layer at the stagnation line, including radiative trans- 
port, was presented by Rigdon, Dirling, and Thomas (1968). This analysis 
included massive blowing (up to 10 percent) for air-to-air injection. In 
1969 they extended the solution to massive blowing of ablation products 
(Rigdon, Dirling, and Thomas, 1970 )> with blowing rates up to 20 percent. 

The numerical procedure which they employed is computationally time con- 
suming. They used an initial value forward integration scheme in which 
they were required to adjust four initial conditions (temperature, tan- 
gential velocity gradient, temperature gradient, and the shear)., 
evaluated at the stagnation point. These initial conditions were iterated 
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until the two boundary conditions for both momentum and energy equations 
were satisfied at the shock and at the body. In the ablation analysis 
they were further required to satisfy the binary species continuity equa- 
tion over the shock layer. If the same initial value forward integration 
scheme is used., then this requires a guess for the elemental mass fractions 
of the injected foreign material and its gradient at the stagnation point 
in order to satisfy the boundary conditions at the shock wave and the 
body. The net result is the requirement that one guess six coupled initial 
conditions which are unknown apriori in order to satisfy the boundary 
conditions. 

From the results of their solution to the exact governing equations, 
Rigdon et al. (1970 ) were able to make direct comparisons with the integral 
results of Wilson and Hoshizaki (l 969) and concluded that the polynomial 
approximations which had worked so well for nonblowing were not sufficient 
to describe the momentum equation solution in the presence of massive 
blowing. Although differences existed in the radiation models (Rigdon et al . 
( 1970 ) used the SPECS code (see Thomas), whereas Wilson and Hoshizaki (1969) 
used an updated version of RATRAP (see Wilson, 1967)), the factors of two 
to four differences in the ablation layer thicknesses could not be explained 
on the basis of radiation model differences. 

In 1969 Wilson ( 1970 ) concluded that the approximate integral solution 
to the momentum equation was inadequate for large mass injection and/or 
the Reynolds numbers of interest. As mentioned previously, his treatment 
of the energy and the species continuity equations was similar to that of 
Howe and Viegas. He used Dorodnitsyn type of transformations (integral of 
the density over physical distance) and similarity conditions to obtain 
exact analytic solutions to the energy and species continuity equations. 
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He attributes the momentum equation solution (see Wilson, 1969) a 
solution technique developed by Chou and Blake for a similar problem. 

Upon performing an additional coordinate transformation which involves 
viscosity and assuming that the density viscosity product is a constant 
across the shock layer, Wilson (1970) developed the second-order differen- 
tial equation. He subsequently differentiated the equation, which makes 
it linear (in the second derivative), to obtain an exact analytical solu- 
tion to the equation in terms of exponential integrals. Rather than 
solve the momentum equation with an initial value technique as was done 
by Howe and Viegas, Wilson employed a successive approximation algorithm 
to all the analytical governing equations. In the successive approxima- 
tion scheme, which is similar to the technique applied in the present 
study, the distribution of properties are initially specified across the 
shock layer and the governing equations are iterated until satisfactory 
convergence is obtained. 

Wilson (1970) observed that with his formulated equations he was 
unable to obtain a numerical solution to the fully viscous equations for 
relatively large blowing rates (greater than 5 percent, approximately). 

As previously mentioned, Wilson traced the problem to one on numerical 
instability associated with taking the differences between exponentially 
large numbers (greater than about e 10 ), which were about the same order 
of magnitude. Since the computer only carries about 8 to 16 significant 
digits, then the resulting difference between these large numbers becomes 
meaningless. Apparently, these exponentially large numbers occur in the 
numerical solution in regions where the viscous effects become small 
(Wilson observed the effect in the inner region near the wall for large 
blowing) and part of the problem could be due to the loss in precision 
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when one uses something like a Simpson's rule (Conte) to numerically 
integrate under an exponential curve. 

Wilson (1970) was able to circumvent the numerical precision problem 
for massive blowing by solving inviscid equations in the inner region 
and the f ull y viscous equations at a somewhat arbitrary distance from the 
body. The interface criteria was that the power to which e is raised 
be less than 10. 

In 1970> Smith, Cuttles, Sullivan, and Graves presented a combined 
flow field and ablation study of a blunt body entering the earth's atmos- 
phere at interplanetary velocities. The analysis, which was motivated by 
a study of a flight experiment to examine the radiation and materials 
response problems at hyperbolic entry velocities, yielded transient 
ablator mass loss rate histories for a complete entry trajectory. The 
results indicated that ablation rates from high-density phenolic nylon 
reached a peak of about O.O55 g/cm 2 -sec (10 percent of free-stream mass 
flow rate) at the stagnation point of a 122-cm-diameter ellipsoid body of 
4:1 axis ratio. Smith, Suttles et al. examined the entire subsonic flow 
region surrounding the blunt body by dividing it into three interacting 
regions: an inviscid outer layer, a boundary layer, and a charring heat 
shield. The inviscid outer layer flow was determined by a one-strip inte- 
gral method with radiation developed by Suttles (1969)- The inviscid flow 
field and the ablation solution provided the boundary conditions for the 
radiating boundary-layer computations. The boundary layer and the 
ablation calculations were iterated until the heating rates and the 
ablation rates converged. For the larger blowing rates the numerical 
method for the solution of the boundary-layer equations was not suitable 
and they adapted an integral procedure by Libby (1962) for the radiating 



l8 


boundary-layer solution. The boundary layer and the inviscid outer layer 
solutions were joined by assuming a cubic variation for the elemental mass 
fraction distribution within the boundary layer and by adjusting the edge 
enthalpy condition. 

One of the significant conclusions of Smith, Suttles at al. was that 
the introduction of ablation products into the boundary layer did little 
to attenuate the radiative flux to the wall. In their analysis they used 
the RATRAP computer code developed by Wilson (1967 ), and they compared 
their results with Chin and with Rigdon et ad. (1969); who used different 
radiative transport computer codes. Whereas Chin's solutions for a 
7.6 percent blowing rate and Rigdon et ad. (1969) solutions for a 
20 percent blowing rate indicated reductions in the wall heat flux of 
about 45 percent below the nonblowing rates. Smith, Suttles ert ad. calcu- 
lated reductions of only about 22 percent. Part of this difference has 
been traced, particularly in the case of Chin's results to differences in 
the radiative transport models. However, Wilson and Hoshizaki (1969) 
in their approximate integral approach were indicating radiative heat 
flux reductions of 40 percent from nonblowing heat fluxes for 10 percent 
blowing rates and of 60 percent for 20 percent blowing rates using the 
RATRAP radiation code. Wilson's more recent results (Wilson, 1970 ), using 
his improved momentum equation solution, have indicated much lower heat 
flux reductions (only 18 percent) at blowing rates of 5 and 10 percent. 
While there are differences in the free-stream conditions in the Wilson 
and Hoshizaki (1969) and the Wilson (1970 ) studies, they do not appear 
to be sufficient to account for the differences in the radiation blockage 
effects in the two analyses. Apparently, the answer to these differences 
must reside in part in the analytical and numerical treatments of the 
governing flow equations. 
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One of the questions which the present analysis will attempt to 
answer is whether the approximate integral treatments of the governing 
equations and/or "exact" numerical treatments of approximate systems of 
equations can accurately define the flow properties required for the 
radiation computations with and without blowing. The numerical technique 
to be developed for the analysis of the coupled viscous, radiating flow 
along the stagnation line and including massive blowing is an implicit 
finite difference scheme. Blottner has shown this approach to be compu- 
tationally superior to initial value schemes for chemical nonequilibrium 
stagnation line studies without blowing. 

In the implicit method, the problem is treated as a two-point 
boundary value problem in which boundary conditions are specified at the 
shock and at the body. The entire shock layer is treated as viscous, 
which requires no "patching" of two or more solutions. The thin shock 
layer equations which govern the viscous along the stagnation line (Ho 
and Probstein) and which are exact through second order are solved at 
discrete nodal points along the stagnation line by iteration through the 
application of a method of successive approximations. Singularities do 
not appear in the formulation since the viscous term takes effect as the 
convective terms (mass flux) approach zero in the vicinity of the 
stagnation point. 
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ANALYSIS 


Thin Shock Layer Equations 

The governing equations for the steady-state flow of a viscous radia- 
ting gas in the stagnation region of a blunt axisymmetric body at moderate- 
to-high Reynolds numbers are given by Ho and Probstein and Scala. They 
are: 

Continuity: 


— — (r'p'u' ) + — — (K'r'p'v 1 ) = 0 
dx' Sy’ 


( 1 ) 


X -momentum: 
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Y -momentum: 
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Energy: 
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(i = 1,...,N species) 

Primed symbols are used to denote dimensional quantities, unprimed 
denote dimensionless quantities. 
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In the body-oriented coordinate system shown la / Figure 3/ the quantity 
, which is the coordinate stretching function, is defined" by 


k* = 1 + K'y ' 


where 


K‘ = K' (x) = -^ 7 - 
and r satisfies the equation 


R b 


(6a) 


(6b) 


dr' = k ' sin 0^ dx' + cos 0-^ dy' 


(6c) 


The quantities q^ , q^ y , and q^ y are the heat fluxes in the 

y' -direction due to conduction, diffusion and radiation, respectively. 

J! is the mass diffusion flux in the y' -direction and cd! is the 

1 

net rate of production of the ith chemical species. 

In comparison to the heat and mass diffusion fluxes in the y'- 
direction, the corresponding fluxes in the x' -direction are generally 
considered negligible (see Ho and Probstein) in the stagnation region. 
These fluxes in the x' -direction are assumed to.be negligible in the 
present analysis also. The thin shock layer equations are the simpli- 
fied boundary -layer equations (Navier-Stokes ) including the curvature 
terms for the stagnation region. The equations are considered to be 


i-l 


(see Ho and Probstein), where 


accurate to the order — r Rey 

P’UIR' \ P s / 

Hey = — when radiation is not considered. 

s u* 

^s 

Stagnation Streamline Equations 

At the stagnation line (x = 0), the conservation equations for the 


thin shock layer reduce to ordinary differential equations upon expanding 
the flow variables in the following power series (based on symmetry 
considerations ) 
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u = a(y)x + a 1 (y)x5 + ... ^ 

v = v(y) + v 1 (y)x2 + ... 
p = P(y ) + P-^y )x 2 + . • • 

H = H(y ) + H-j_(y )x 2 + ... 
a ± = a ± (y) + a i;L (y)x 2 + ... 


( 7 ) 


and applying the geometry relations (eq. (6)), then taking the limit as 
x -» 0. The limiting forms of the governing equations for the stagnation 
line flow become: 

Continuity: 
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X -momentum: 


Y -momentum: 
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Energy: 
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Species continuity: 
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where 
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( 13 ) 


P = 3(y) = - 


>&x 2 /x=0 


( 14 ) 


The foregoing conservation equations have been written in nondimen- 
sional form from the following set of dimensionless quantities: 
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where the subscripts 00 and s refer to free-stream and post-shock 

conditions, respectively. 

Restrictions and Assumptions 

The governing stagnation line equations are general equations and 
are restricted only by the requirements that P g /Rey s « 1 and that the 
radiative heat flux in the x-direction are comparatively negligible. At 
this point it is desirable to establish the state of the gas and the 
diffusion and heat flux models to be used in the study. The following 
basic restrictions and assumptions are imposed on the governing equations: 
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(1) The gas is in local thermodynamic and chemical equilibrium. 

( 2 ) Diffusion of the chemical species is governed by a binary 
diffusion model. 

( 3 ) Radiative energy transport occurs within a one-dimensional, 
infinite planar slab (tangent slab). 

Anderson indicates velocity and altitude limits for equilibrium and 
nonequilibrium chemistry, and it is observed that for the conditions of 
interest in this study, chemical equilibrium is considered to exist. 
However, to date no detailed calculations which consider the ablation 
products have been made to firmly establish this assumption. 

The binary diffusion model, which assumes that the individual chemi- 
cal species have the same mobility as the two chemical species which are 
used to represent the binary diffusion process, has been assumed in all 
previous radiation studies of this problem. 

For the blowing situation, it is assumed that all the numerous 
chemical species of the ablation products can be represented by a single 
chemical specie and all the chemical species of air can be represented 
by the other single chemical specie to form the binary diffusion coef- 
ficient. Rigdon et al. ( 1969 ) examined atomic carbon-atomic nitrogen 
diffusion and atomic hydrogen-atomic nitrogen diffusion and concluded 
that the actual choices of the two chemical species to form the binary 
diffusion coefficient did not appear to be critical in the solution. A 
similar analysis of the effect of diffusion coefficient on the radiative 
heat fluxes will be made in this study. 

Wilson (1970) indicated that he was considering a pseudo multi- 
component diffusion analysis (bifurcation model, see Graves) in which 
the individual binary diffusion coefficients computed for each combina- 
tion of two individual chemical species are statistically lumped to 
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form an effective binary diffusion coefficient; however, the calculations 
have not been made to date. 

The errors expected due to the tangent slab approximation have been 
previously discussed in the review of the literature and are about 
5 percent or less for this problem (see Kennet and Strack, Koh and 
Hoshizaki, and Lasher). This approximation is made in RATRAP for the 
computation of the heat flux term that appears in the governing energy 
equation (eq. (ll)). 

Although additional assumptions, which are subsequently discussed, 
are required with regard to boundary conditions and second-order transport 
effects, the foregoing three assumptions are the most basic and restrictive 
assumptions which are required in this analysis. There is no way to estab- 
lish firmly the validity of these basic assumptions for this problem, beyond 
the studies which have been previously cited, without introducing addi- 
tional complexity to the calculations . However, it is anticipated that 
the implicit finite difference scheme should lend itself readily (con- 
ceptually) to both chemical nonequilibrium studies as demonstrated by 
Blottner's (1970 ) analysis and to multicomponent diffusion studies (Graves). 
It should be particularly useful for the nonequilibrium studies because 
stability requirements are not as stringent in this approach as in explicit 
forward integration schemes. 

In order to uncouple the stagnation line solution from the remaining 
subsonic field, it is necessary to assume the relationship between the 
shock and the body curvatures at x = 0. In this analysis, as in most of 
the prevalent analyses, it is assumed that the shock and body are concen- 
tric. The results of Suttles (1969) inviscid radiating analysis indicates 
that the assumption is reasonable. 
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Preheating of the ambient air upstream of the shock wave due to 
radiation transport (precursor effects) is neglected. Smith, Hoshizaki, 
and Lasher, and Rigdon et al. indicate that precursor effects begin to 
become significant at velocities around 17 to 18 km/sec and above. Thus, 
the present solutions, which employ the Rankine-Hugonoit conditions (see 
Hayes and Probstein) for the discontinuous jump conditions across the 
shock, will be restricted to velocities somewhat lower than this. 

In the numerical solutions, unless otherwise specified, a Newtonian 
pressure distribution (see Hayes and Probstein) is used to evaluate (3, 


the 


\ctx 2 / x=0 


term in the x-momentum equation, i.e., 


8x 2 / v=l 


= - 2 . 0 . 


x=0 


It should be noted that Wilson (1970 ) used a value of -3-0 which can lead 
to a thinner shock layer in his calculations. This effect will be inves- 
tigated in the present study. 

The radiation transport computer code which is used in the radiation 
computations is RATRAP, developed by Hoshizaki and Wilson (see Wilson, 1967 ), 
which considers most of the primary radiating chemical species associated 
with carbon, oxygen, nitrogen, and hydrogen mixtures. The detailed thermo- 
dynamic and chemical composition calculations for equilibrium chemistry 
are performed in the computer code FEMP, developed by Brown et al. FEMP 
is included by Wilson as an integral part of RATRAP. 

In the analysis, it is assumed that the transport properties for air 
developed by Hansen apply for both air-to-air injection and ablation - 
products-to-air injection. This assumption for the ablation products will 
be superficially analyzed by perturbations in the pertinent transport prop- 
erty (Prandtl number) to determine its effect. It should be noted that 
Rigdon et al. ( 1969 ) ran two identical cases with the exception of pure 



air transport properties in one and the combined ablation products air 
transport properties in the other , and the resulting differences in the 
radiative heat fluxes to the wall were less than 1 percent. 

The continuity and y -momentum equations (eqs. (8) and (10)) are 


unchanged when the foregoing assumptions are applied. The only perturbation 

'82p\ 

will be 


in the x -momentum equation (eq. (9)) is the j- term arising from the 


0 


Newtonian (or specified) pressure distribution in which 
taken to be a constant in y in this analysis. 

The heat flux terms in the energy equation (eq. (ll)) can now be 
defined by applying the binary diffusion and tangent slab approximations. 
The conductive heat flux term is given by Fourier's law of heat conduction 


, , , dT’ 

q c,y = - k <= 


(16) 


The multicomponent diffusion heat flux is given by Bird, Stewart, and 
Ligjitfoot as 
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where l — 1, ... ,N species . 
For binary diffusion, 


da,- 

j; = -P'D' —4 

1 12 dy' 


(18) 


For constant pressure (pressure varies about 2 percent across the 


shock layer for this study ), 
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Equations (17), (18), and (19 ) can be combined to yield 



The term in brackets of equation (20) has been defined by Hansen to 
be the "reactive" conductivity , k-y which yields a form of Fourier s law 
of heat conduction that describes the energy transport by binary diffusion. 
Hansen has computed and tabulated kp and the ’ lumped or effective 
conductivity k = k^ + kp for air system for temperatures up to 15,000 K 
and pressures up to 100 atmospheres. 


The radiative heat flux equation for the heat flux at a point within 
a one -dimensional slab is developed in the Appendix. That equation is 



and is described by Kourganoff. B y is the Planck function and a v is 
the modified linear absorption coefficient which is a function of the 
temperature, pressure, and the chemical species number densities within 
the slab. The equation is valid for a non-gray self -absorbing gas. The 
absorption coefficients and radiative heat fluxes are computed in the 
radiation computer code RATRAP developed by Hoshizaki and Wilson (see 
Wilson, 1967). The radiation model takes into account both continuum and 
atomic line radiation exchange in a slab of non-uniform temperature. 
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The simplified energy equation can now be written as 
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( 22 ) 


where q' , is given by equation (21). The combined conduction- 
diffusion heat flux term can be written in terms of total enthalpy from 
the following relations for static, total enthalpy, and "effective" 
Prandtl number: 
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Upon applying the nondimen sional relations defined in equation (l4 ) and 
rearranging equation (25), the dimensionless energy equation becomes 
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The species continuity equation (12) undergoes considerable modifi- 
cation because of the binaiy diffusion and equilibrium chemistry assump- 
tions. The species continuity equation was given by 


0 dix-? 
K^PV — — 
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Under the local chemical equilibrium assumption, the volumetric rate 
of production of species i is indeterminate. However, the fact that 
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the net rate of production of the mass fractions of the chemical elements 
is zero can be utilized to yield a tractable solution to the species con- 
tinuity equation. Assuming that the binaiy diffusion law holds for the 
mass fractions, jL.£. , 


_ do,.- 

J i = - pD i2 17 


( 27 ) 


where the barred quantities refer to elemental mass fractions, the species 
continuity and diffusion equations can be combined to yield the elemental 
diffusion equation 
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The dimensional binary diffusion coefficient is given by Hirschfelder 
et al . as 
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where 
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c • 

The molecular constants a± and — i are tabulated by Svehla and are 

K 

given in Table 1 for the chemical species of interest. The reduced 
collision integral, is based on the Lennard- Jones (12-6) poten 

tial and is taken from Hirschfelder et al . as a function of the nondimen- 
sional T^ -which is given by 



For the binary diffusion model assumption, the individual elements 
of the ablation products and of the air which passes through the shock are 
considered to diffuse in the same respective manner as the two chemical 
species (one for the ablation products, the other for the air) used to 
form the binary diffusion coefficient. Thus, a distinction need not be 
made between the individual elements but only between the total mass 
fractions which represent the ablation products and the remaining mass 
fraction which represents the air products. Since the total mass fraction 
of all the elements (and, for that matter, the chemical species) must 
equal unity at any point, then one need solve only one elemental diffusion 
equation which is given by 


_d_ 

dy 


K 2 pDi 2 



o dotp 
hc^pv — — = 0 

dy 


(30) 


where dp 
products, 
puted by 


is the total mass fraction of the elements of the ablation 
The total mass fraction of the air products, a A , is then com' 






( 31 ) 
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The mass fractions of the individual elements are then calculated by the 
equations 


£ J-j rnj 


i = 1,N 


( 32 ) 


where a. . . is the specified elemental mass fraction of the ith ele- 

lj, mj 

ment of the ablation products that is injected at the wall and M is 
the same ith element that basses through the shock layer. The density 
and the individual chemical species are then calculated from an equation 
of state by the energy minimization subroutine (FEMP) where 


p = p(a.,h,p) 

( 33 ) 

a i = a i (aj,h,p) 

For viscous radiating stagnation line analyses, the resulting 
governing equations (8), (9), (10)., (26), and (28) are the most general 
system of equations that are treated in the literature. These are the 
exact thin shock layer equations (for the tangent slab and binary diffusion 
assumptions) for the stagnation line which are solved in this analysis. 

As previously mentioned, Rigdon et al. (1968) and (1969) have applied 
initial value techniques for the numerical solution to this system of 
equations (with the exceptions that they assume constant pressure across 
the shock layer and neglect the curvature terms) for the massive blowing 
problem. 

Transformed Equations 

Before proceeding to the development of the solution technique, it 
is desirable to transform the governing equations by a change in independ- 
ent variable from y to q, where 




( 34 ) 
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This transformation has the important effect of fixing the shock boundary 
at rj = 1. The inclusion of density in the transformation gives, for a 
fixed nodal spacing in t), a finer resolution on a physical scale, y, 
at points near the body and within the boundary layer where conditions 
are rapidly changing. 

The transformed system of equations become: 

Continuity: 


d(K 2 py) 

dtj 


25 Ka 


(35) 


X -momentum: 


jiaL kA 

Rey s S li - 'l ' iy 


2 

Kp V da 
5 dr) 


pa 2 + 


Kpva 


-3 


(36) 


Y -momentum: 


dp dv 

— = -pv — 

dTj dq 


(37) 


Energy: 



fr 2 pp 

dr)\ 

, Pr 


dr) _ 


? dH 

+ 5Rey s ^pv — = 


dA 2 | 

dT)\^ Pr 


^ Pv S ) + 5 Re Ys^( K2 l 


dT), 


'dn 


■R,T) 


) (38) 


Elemental diffusion: 


_d_ 

dr) 


k 2 P 2 D. 


12 



(39) 


Equation of state: 


p = p( a;i ,h,p) 


fro) 


The reader is cautioned that the normal velocity, v, has been 
redefined as negative in the positive y or r) direction in the 
above equations. This change in no way affects the solution to the 



equations since the boundary conditions, as will subsequently be discussed, 
have been appropriately modified to reflect this change. The primary 
reason for the change in the sign of v is to aid future users of the 
computer program to understand the signs on the equations as they are 
programed. In programing the governing equations, the author chose to 
define velocities in the direction of the free-stream flow (negative 
y-direction) as positive. The flow field coordinate system for the stag- 
nation line in the transformed coordinates is shown in Figure 4. 

Eta is specified as zero at the body and as unity at the shock. The 
normalizing parameter, 5, is an unknown in the problem which is obtained 
from the solution of the continuity equation (35 ) and is given by the 
relation 


6 


-(k 2 Pv) w + (k 2 Pv) s 



(41) 


Boundary Conditions 

The subsonic flow field in the nose region is governed by elliptic 

equations; consequently, the stagnation streamline solution is influenced 

by the flow within the entire subsonic region. However, this influence 

only enters in the (3 term, i.e. , (—£) in the x-momentum equa- 

tion ( 3 6 ) and in the curvature of the shock wave which provides the 

boundary condition for the tangential velocity gradient, as at the shock 

wave. For a shock wave and body which are concentric, the Rankin-Hugonoit 

/d 2 p\ 

relations yield the following relation between — and the velocity 

\ cbt / x=0 

gradient, a s (see Hayes and Probstein) 



35 


P = 


VSx 2 / x=0 





Newtonian impact theory predicts (see Hayes and Probstein) 


Sx2y x= o 


- 2,0 


(42) 


(43) 


thus. 


ifpy /2 = x.o 

2 c 


With 3 and a s specified, the stagnation streamline problem 
becomes a two-point boundary value problem, where conditions are specified 
at the shock wave and at the body. The boundary conditions necessary to 
solve the transformed conservation equations for the flow of a viscous 
radiating gas with foreign species injection, where the chemical species 
are in local thermodynamic equilibrium are: 

At the body surface (wall), rj = 0 (y = 0): 


a = aw = 0 

pv = (Pv) w 

P = P w « P s + -(PV ) s 

H = Hw 

a. = a. . . (and 54, = 1 ) 
x i,mj * 




) (44) 


i = 1, ...,N elements J 


where the' 5^=1 boundary condition applies for condition of strong 


blowing. 
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At the shock wave, H = 1 (y = y s ) : 

B* = Q*g = 1.0 

pv = (pv) s = 1.0 

P = P s 
H = H g = 

a- = a ±}00 (and 0 ^ = 0) i = 1, ...,M elements 

where a,- is the mass fraction of element i which is injected at 

r, xnj 

the wall, and cLj_ «, is the mass fraction of element i which passes 
through the shock wave from the free stream. Conditions immediately 
behind the normal shock are computed from Rankine-Hugonoit relations. 

Heat fluxes across the slab boundaries are assumed to be zero, i. e. , there 
is no precursor heating of the free stream by the shock layer and no radia- 
tion from the body into the shock layer. 

Numerical Solution 

The finite difference approximations to the governing equations 
(35) to (39) and the procedure for the numerical solution to this sys- 
tem of equations is developed in this section. The equations are 
written in finite difference form for a network of N equally spaced 
(in q) nodal points between the body (n = 1) and the shock (n = N) 
which are shown in Figure 5 . For most of the calculations, 21 nodal 
points are used. The overall numerical solution technique is iteration. 
Profiles and parameters are assumed initially for each nodal point 
across the shock layer. The governing finite difference equations are 
then solved sequentially, using the most recent values of the profile 
parameters until satisfactory convergence is obtained at each nodal 



point. 



37 


Finite Difference Approximations 

Either two- or three-point finite difference formulas are used to 
numerically differentiate and integrate the governing differential equa- 
tions. Where possible, three-point central differences are employed 
since they are accurate (at a particular nodal point) to order (Aq) 
where At] is the distance between points, whereas two-point differences 
are accurate only to order Atj (see, _e.g_. , Conte or Crandall). Two- 
point windward differences (with the flow) are employed for the convec- 
tive terms where dictated by stability requirements of particular 
governing equations. The stability requirements are discussed in the 
Results and Discussion chapter. 

The tabulation of the finite difference formulas that are employed 
to approximate the derivatives is given below (see, £.£• , Conte or 
Crandall for the development). The formulas are valid for equally 
spaced increments in At). 

The central differences formulas are: 


df \ = - f n-l + f n+l 

\ dq/ n " 2 AT) 


(46a) 



f n-l ~ 2f n + f n+l 
(Ati) 2 


(46b) 



" G n^n-1 ” ( c n-l ” c n+l^n + c n^n+l 
2 Ati 


(46c) 


A 

dri 



( c n-l + c n^n-l ~ ( c n-l + ^ G n + c n+l^n + ( G n + c n+l^n-il 

2 (Ati ) 2 


(46d) 


where c and f are arbitrary functions. 



In developing equation (46d), it has been assumed that c varies 


linearly between nodal points. 

The windward difference formulas are either forward or backward 
differences^ depending on the direction of the flow. If the mass flux, 
pVj is positive (toward the body); forward differences are taken for the 
convective terms; and if pv is negative^ backward differences are taken. 
The formulas are: 


If (Pv) n > 0, 


If (pv) < 0; 
n 



(47) 


(48) 


Numerical integration is performed by the Simpson's rule approxi- 
mation (Crandall) over parabolic sections of the profiles. The finite 
differences equations used in performing the quadrature over the interval 

from Ti to T) , and over the full interval from 1} to ri „ are: 
n n+± n u *- 


h+1 ^ At) 

f (t)^ 7 ! = — 


•n 


5f n + 8f R+ 2 “ I n +2 




and 


(49) 


\+2 . _ AT) 

f(ri)dT) = — 

i 3 


f + 4f + f , 0 

n n+l n+2 


'n 



Initial Profiles 


For the general case, the following profiles are initially assumed: 


Pv(ri) = (pv) w + 


(PV) S - (pv)„ 


(50a) 


p(n) = p s + |(pv 2 ) s ( 1 - T i 2 ) 


(50b) 


a(t] ) = a s rj 


(50c) 


h ( r] ) = h w + (h g - h w )r] 


(50d) 


^(n) = (£p ) w (! - T i) 


(50e) 


p(r]) = p(a i (n),p(r,),h(ri)) (50f) 

The density is computed in FEMP for a gas mixture in chemical equilibrium. 
All other quantities appearing on the right-hand side of the initial pro- 
file equations are available either from the input boundary conditions or 

they are computed from Rankine-Hugonoit relations for the normal shock. 

P 7 ! dh 

The scale factor K=l+Ky=l+K 6/ — which, for large 

J 0 P 

Reynolds numbers, typically is taken to be a constant of unity across 
the shock layer (Ho and Probstein) has been retained in the governing 
equations for this analysis. While k exhibits only weak variations 
across the shock layer (for the typical cases of interest in this analy- 
sis k varies from about 1.0 to 1.02 from the body to the shock for no 
blowing and from about 1.0 to 1.1 from the body to the shock for strong 
blowing (20 percent of (pv) ), the computational time required to 
generate k is insignificant in this problem. However, an initial 



value for the transformed shock displacement, 5, is required as a 
consequence of retaining the k. The physical shock displacement dis- 
tance is correlated hy Inouye to "be 


for no blowing. 

density of p* 
s 


Since 


y 


R b 


0.78^ 

R b p s 


5P “J 


r 


— , then for 
P' 


a constant shock layer 


6 ~ O.78 


For moderate blowing, 

S ~ 1 (51) 

which is the initial value assumed in this analysis. 

This completes the statement of the problem of determining the 
initial profiles for the general case. The appropriate boundary condi- 
tions and initial profiles for special cases of interest, such as no 
blowing or air-to-air injection, are covered in the subsequent dis- 
cussion of these special cases. 

Solution Procedure and Finite Difference Equations 

With the assumed initial profiles and parameters (eqs. (50) and (51)), 
the conservation equations are solved by successive iteration. The 
coupled system of equations are solved in the following sequence: 
continuity, y-momentum, elemental continuity, x-momentum, energy, and 
equation of state. The most recent values of the profile parameters 
and 6 are used in the computations. The density distribution is 
necessary to solve the governing differential equations. After these 
equations have been solved, updated values of density are computed from 
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the equation of state which is compared with the previous density values 
at each nodal point. The entire sequence through the governing equations 
is repeated with the new density distribution until two successive passes 
yield nearly identical (within 1 or 2 percent) density values at corre- 
sponding nodal points. The flow diagram which illustrates the solution 
procedure is shown in Figure 6. 

With the distribution for the tangential velocity gradient and the 
scale factor k(t]), the continuity equation (35) is numerically inte- 
grated by Simpson's rule (eq. (49)) to give the updated shock displace- 
ment distance, 8: 


8 


.(pv) w + (k ov) s 


Ka dTj 


0 


and the mass flux distribution Pv(tj): 


(52) 




Pv(t]) = 25 j *a + (Pv) 
J 


(53) 


The velocity profile, which is required in the solution of the 
y-momentum equation is computed from the previous density distribution, 
p (rj ) , and the updated mass flux distribution, pv(rj), namely 


v(q) = 


Pv( T l) 

P(4) 


(54) 


The y-momentum equation (33) is numerically integrated by Simpson's 
rule (eq. (45)) to give the pressure distribution p ( rj ) : 
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P 



dv 

pv — dT) + 

dr] 


P S 


(55) 


The Pv(r|) and v(r)) profiles are tabulated from the continuity 

equation solution and a central difference numerical differentiation 

scheme (eq. (46a)) gives — ^ ^ — — at the N— 2 points within the shock 

laver. The values of — — at the body and the shock boundaries are 

J dy 

obtained from two-point forward and backward difference approximations, 
respectively. 

The elemental diffusion equation (39) is cast in finite difference 
form by using the central difference formula (eq. (46d)) for the second 
derivative diffusion term and the windward difference equations ( ( 47 ) 
and (48)) for the first derivative convective term. The resulting equa- 
tions are: 

For n = 1, 

(S^O-l = 1 (56a) 


For n = 2, . • - ,N-1, 


(K 2 p 2 D 12 ) n _;L + (xVD 12 ) n - 25 (AT) ) ( K^pv ) n (l - SIGN) 


2 2 . 
‘ *— 


(a F ) n _] 




+ 25(AT])(K 2 pv) n (l - SIGN) 


(SfO n + 


(K 2 p 2 D 12 ) n + (^^ 12^+1 


2 ^ 2 t 


+ 26 (At) ) ( ^ Pv ) n ( SIGN ) 


^n+1 0 


( 56 b) 



For n = N, 
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(£p) N =■ o (56c) 

where 

SIGN = 0 if (pv) n < 0 

and 

SIGN = 1 if (pv) n > 0 

The formulation results in a handed (tridiagonal) matrix system of 
equations of the form 



where B n , and C n are the coefficients of the (a F ) n _q> (a F ) n > 

and (cL ) n+1 terms, respectively, and D n are the terms appearing on 
the right-hand side in equation (5 6 )• At the wall boundary, n = 1, 
the coefficients are Bq = 1, Cq = 0, Dq = 1 and at the shock boundary, 
n = N, the coefficients are Ajq = 0, Bpj = 1, Dim = 0- The system °f 
equations is easily solved by Potters' method, which is a form of 
Gaussian elimination that is efficient for solving a banded matrix system 
of equations (see Potters and Conte ) . The solution to equation (56) 
yields the total mass fraction distribution, at N nodal points, of the 
ablation products, subject to the constraint that 0 < (a F ) n < 1* The 



44 

mass fractions of the individual air and ablation product elements are 

then computed from equations (31) and (32). 

The transport properties , p and P r , required in the solutions to 
the x-momentum and energy equations are obtained from a tabular lookup 
as functions of the temperature, calculated in the FEMP subroutine, and 
the pressure, obtained in the y-momentum equation solution. 

The x-momentum equation (36) is next solved. The equation is cast 
in finite difference form using the central difference expressions given 
by equation ( 46 ) to approximate the first- and second-order derivatives. 
Since the equation is nonlinear in the velocity gradient, a, an iterative 
procedure is required for its solution. A quasi -linearization approach 
(see Bellman and Kalaba) is employed which generates the solution by a 
rapidly converging iteration. 

The x-momentum equation is quasi-linearized by the following tech- 
nique: At a given nodal point, let 



a 


(i-D 


(i) 

- a' 1 ' 1 ’ 

+ 

a (l) - a 


- 


(58) 


where the superscripts i and i-1 refer to the values of a(r}) at 
the ith and i-lth iterations, respectively. 

Upon assuming that 


lim 

i — » 00 


a^-a* 1 - 1 ]' 


( 59 ) 


the x-momentum equation becomes linear in a^(t]) where 


L (i) 


2a 


(i-l)„(i) 


■ M 


as i -» °° 


(60) 



Combining the finite difference formula for the derivatives 
(eq. (46)) and equation (60) yields the finite difference form of the 
x -momentum equation: 

For n = 1, 

aj 1 ^ = 0 (6la) 


For n = 2 , . . . ,N-1, 


( K P 2 v) n _ ( Kp )n[(PP)n-l + (PP)nj \ (i) + 


2S(Ari) 


25Rey s (Aq) £ 


a n-l + < 2 


pa 


(i-D 


(Kpv) + 
v n 


(Kp) n [(pq) n _! + 2 (p|i ) n + (Pd) n +ljl (i) 


25Rey q (Ar]) 2 


n 


(K P 2 v) n (K P ) n Qpn) n + (PH) n+ j\ (± 


128 (Aq) 


28Rey c .(ATj)' : 


(i) 

n+1 


= 0 + <P 


(i-i)l 2 


(6lb ) 


n 


For n = N, 



(6lc) 


The subscript n refers to the nth nodal point in the shock layer and 
the superscript i refers to the ith iterative value of a. 

For a linear iteration on the velocity gradient, a, let 



a i_1 + e 
n n 


(62) 


where e n is the error in the ith iteration at the nth nodal point. 
Substituting equation (62) for a^ into equation (6l) yields a tri- 


diagonal matrix system of equations of the form: 
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Vn-1 + Vn + C n £ n+1 " D n n ' ’ ^ N_1 


( 63 ) 


where the A^, B n , and C n are the coefficients of the a^^, a^' \ 

and a^^ terms, respectively of equation (6lb) for n = 2, ...,N-1. 
n+1 

The term D n is given by 


D n = p + < P 


a 


(i-1) 


(i-1) t> „ Ji-1) 


- Vn-1 - B n a i 


n 


C n a n+1 


( 64 ) 


n 


for n = 2, . ..,N-1. Since the boundary conditions are specified at the 
body and at the shock, then the equations at the boundaries become: 


61 = 0 (or At = Cp = Dp = 0 and Bp = 1) 


e ]\f - 0 (or Ajj - Cpj- - - 0 and 1) 




(65) 


J 


The matrix system of equations (64) and (66) is solved by Potters’ 
method and is iterated until 


e n < T 


¥• n 


In this study, this assigned convergence interval, t , is 0.01. 

The total enthalpy distribution is obtained from the tridiagonal 
matrix solution to the energy equation (38) in a single pass. As pre- 
viously discussed, the total enthalpy derivative for the combined convec- 
tive and diffusive heat flux terms is included explicitly in the matrix 
solution for the H ' s, whereas the radiative flux term is computed using 
the enthalpy profile of the previous iteration and the updated pressure 
and elemental mass fraction profiles from the y-momentum and species 
diffusion equation solutions, respectively. . The three-point central 


difference formulas for the derivatives (eq. (46)) yield the following 
finite difference form of the energy equation: 


For n = 1 , 


H 1 = Hw 


(66a) 


For n = 2., . . . ,N-1* 


1 


2 Arf 


+ [l U 

V p r / n-1 \ P r /n 


~ “(^)n H n-1 

2 AT) j 

-ih 


|2(Atj) £ 


* 2 pA + 2 r 2p A + (— ) 

v 2 r /n-1 \^r Jn \ p r /n+1 


H 


n 


2 Arf 


^ 2 PPj ^ / K 2 p|j. \ 

V P r / n \ P r / n+1 


+ ^(k 2 pv) ) H ^ 
2 At) n+ ^ 


\ 


— — 

d / k 2 li dv 1 

+ 5Rey 

d f o \ 

— — Z pv — 


dT)\P r dT)/ 

n 


- -* 


n 


(66b) 


For n = N, 

% = n» < 66c > 

The derivative terms appearing on the right-hand side of equation (66b) 
are evaluated by the central differences formulas (eq. (46)). 

The density distribution which was required in the solution to the 
governing differential equations can now be updated by the equation of 
state 


P n = P 


n KA' 


P n >hn 


(67) 


at each nodal point. The updated values of P n provide the mechanism 
for iterating the governing equations. Upon setting P n equal to either 



the newly computed values or some percentage of the old values of density, 

and the new values p^^ (in order to speed convergence of the 
n 

flow field solution), the entire procedure, beginning with the continuity 
equation, can be repeated until 

p (i) _ p (i-l) I 

— n _n l_ < e ¥ n ( 58 ) 

p(i) 

*n 

The interval of convergence e used in this analysis is 0.02 or less. 
Salient Features of the Implicit Finite Difference Algorithm 

The implicit finite difference algorithm developed in the previous 
section is really quite simple, yet sufficiently flexible to treat the 
viscous radiating shock layer problem. 

With the two- and three-point difference approximations for equally 
spaced increments, the governing thin shock layer equations can be made 
amenable to numerical solution without making unduly restrictive assump- 
tions for the purpose of yielding analytically tractable solutions. The 
tridiagonal matrix system of equations which results can be efficiently 
solved by Potters' method, which requires only about 3 n computations, as 
contrasted to the (n) computations required for a full matrix 
inversion. 

Simple linear and quadratic profiles can be initially assumed as 
functions of the boundary conditions. Or, one may take advantage of 
prior knowledge of the solution behavior to begin with improved estimates 
of the profile parameters (i.._e., read in the parameter values at each 
nodal point). 

No singularities appear in the governing finite difference equations 
since division by pv is not required in the present formulation of the 



problem. In the region where Pv (or v) approaches zero, the diffusion 
terms become important in the elemental diffusion equation (eq. (5 £))• 
Likewise, the viscous terms predominate in the x-momentum and energy- 

equations (eqs. (6l) and (66)) near Pv - 0. 

The method is an implicit scheme in which the unknown quantities at 
point n are calculated as functions of conditions at surrounding points 
as well as conditions at the point itself. This is in contrast to 
explicit schemes where the unknown quantities at a point are evaluated 
solely as a function of conditions at a former point (such as the func- 
tion and its derivative being evaluated at the n-1 point and then 
extrapolated to the nth point to determine the unknown function). For a 
given step size, implicit schemes are unconditionally stable (bounded), 
whereas explicit schemes may or may not be stable. 

In the present approach to the two-point boundary value problem, 
the boundary conditions are specified at the shock wave and at the body. 
Since the matrix system of equations (for a particular governing equation) 
is completely coupled across the entire shock layer (e.g., note the 
appearance of (3^ in the n-lth, nth, and n+lth elemental diffusion 
equations (eq. ( 5 6 )) and the known boundary conditions are included in 
the system of equations , then the computed distributions between each 
endpoint are bounded and are monotonic in behavior (with the exception 
of the pressure, which reaches a maximum at the stagnation point). Thus, 
one has near-maximum information about the behavior and levels of the 
unknowns which are to be computed and also a high degree of assurance 
that the computer run will not be aborted because of overflow (numbers 
larger than the computer will accept). 



50 


Perhaps one of the biggest advantages of the present approach over 
initial value methods resides in the limited number of unknown quantities 
which must be evaluated in order first to satisfy the governing differ- 
ential equations and, second, to provide the necessary inputs for the 
radiative flux computations. In the present method, the unknowns are 
only the properties themselves (pv, p, etc. ) at each nodal point about 
which one has maximum information as to the bounds on the values and 
the general behavior of the properties across the shock layer. In con- 
trast, consider the computations which must be performed in an initial 
value treatment of the problem (such as Runge-Kutta forward integration). 
In the initial value approach for the solution of a single governing 
equation, the unknown property and also its derivative must be computed. 
Generally, little information is available as to the behavior of the 
derivative across the shock layer. As a consequence, the function or 
derivative changes from point to point are generally closely controlled, 
by restricting the step size, to maintain stability as the calculations 
proceed downstream. 

Yet there are several clear advantages of the initial value treat- 
ment of the complete system of coupled governing equations. If all 
boundary conditions are matched at the downstream side in the iterative 
process, then, because of the strict stability requirements, one is 
reasonably assured that all governing equations have been satisfied 
both locally and globally, and the detailed distributions are sufficient 
for the radiation computations. Whereas, in the implicit finite differ- 
ence approach, there is no assurance apriori that the coupled system of 
equations will converge in the overall iteration scheme. It may be 
possible that two equations (_e.£. , the continuity and the x-momentum) 
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may interact with one another such that the pv and a profiles 
oscillate back and forth and never converge- Further, there is no con- 
trolled step size from stability considerations. Thus, care must be 
exercised in inputting the nodal spacing which will satisfy the govern- 
ing equations both locally and globally. An indication of the required 
nodal spacing can be obtained by a cursory examination of stability 
requirements for the individual equations and by decreasing the step size 
until the numerical solutions become asymptotic. These aspects, along 
with the convergence behavior of the overall solution, are examined in 
the Results and Discussion chapter which follows. 
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RESULTS AM) DISCUSSION 

In this chapter the numerical results obtained from the solution to 
the governing equations developed in the previous chapter are presented 
and discussed. The solutions which were obtained in the various phases 
in the development of the method are presented in the sequence of their 
development, beginning with the constant density solutions and continu- 
ing through the viscous radiating shock layer solutions with ablation 
products injection. Resuits are compared with results from existing 
approaches to the radiating flow field problem. 

All numerical solutions were generated on a Control Data Corporation 
CDC 6600 digital computer. 

Constant Density Solutions 

The results of the constant density study are presented in this 
section. The convergence behavior of the solution for the reduced 
system of equations that govern the flow of a constant density gas is 
examined. Selected solutions are shown for inviscid and viscous flows. 
Results are presented for a range of blowing rates and Reynolds numbers. 

The overall solution to the governing equations is by iteration 
until satisfactory convergence is obtained on the density. For a 
constant density assumption (and a viscosity and Reynolds number speci- 
fication required for the viscous flow solution), the continuity, 
x-momentum, and y-momentum equations can be solved independent of the 
energy and diffusion equations. Naturally, the density is not iterated; 
however, the remaining system of equations (continuity and x- and 
y-momentum) are still coupled and, in the present approach, an iterative 
scheme is required for their solution as indicated by Figure 6. The 



equations are solved "one at a time" by a method successive approxima- 
tions until the computed values of pv, a ., and p for two consecutive 
iterations are within a specified accuracy at each nodal point. The 
constant density solutions thus serve a twofold purpose. First , the 
solutions provide a most fundamental means of studying the convergence 
behavior of the coupled system of equations for a wide variety of blow- 
ing rates and Reynolds numbers and thus permits an assessment of the 
adequacy of the successive approximation iterative scheme. Second, it 
provides an indication as to the nodal spacing requirements for the 
implicit finite difference equations formulation. 

Several cases were run for a range of typical blowing rates 
((pv). r = 0 to -0.2), Reynolds numbers (Rey c = 1 to lCp), and nodal 
spacings (N = 11 to 101 ) for a constant density of 20. and a viscosity 
of 1.0. All solutions required about three or four iterations and a 
total 2 to 4 seconds of computer time. No discernible differences were 
observed in the computed pv, p, and a profiles or 5 for N = 11 
and N = 101. It was observed that the computer times were a function 
on the number of nodal points (the 2- and 4-second run times were for 
the N = 11 and N = 101 cases, respectively), whereas increasing the 
blowing rates or Reynolds number had little or no effect on the com- 
puter run times. 

The convergence behavior of a typical constant density is shown in 
figure 7 where the successive solutions to the continuity and x-mamentum 
equations are plotted. The solution is for a constant density of 20. , 
a constant viscosity of 1.0, a typical Reynolds number of 10^, and a 
blowing rate at the wall of -0.2. Twenty-one nodal points were used in 
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the computations. Initial linear profiles for Pv(tj) were assumed. 
The accuracy criteria required that the computed values of pv(n); 
p(r)), and a(r|) for two consecutive iterations be within 0.01 at each 


nodal point (i.e.; 


PV 


(i) 


Pv) 


(i-1) 


<0.01 ¥ n). This accuracy 


n n 

criterion was maintained throughout the entire course of the present 
study of viscous radiating flows. As indicated by the legend in 
Figure 7 , the solution converged in four iterations which required 


about 3 seconds of computer time. 

Three constant density solutions for no blowing are shown in 
Figure 8: an inviscid case and two viscous cases for Reynolds numbers 

of 1 and 1(P. The Rey s = 10^ case is typical of the Reynolds number 
of interest^ whereas the Rey s =1 is actually outside the limits of 
applicability of the thin shock layer equations and is shown merely to 
demonstrate that the present approach yields solutions over the entire 
range of Reynolds numbers. The effect of the boundary layer and its 
extent can be seen in Figure 8(a). For inviscid flow there is no 
boundary layer and a nonzero velocity gradient exists at the wall which; 
from equation ( 3 6 ); is given by a-j_ = J p/Pp- For the Rey s = 1°^ 
case, most of the shock layer is inviscid; with only a thin boundary 
layer present near the wall (to tj » 0.03) where the velocity gradient 
slope in y (i.e., da/dy) is a maximum; whereas for the Rey s = 1 
case, the entire shock layer is viscous and no abrupt changes in da/dy 
are observed. The plot in Figure 8(b) indicates that there is little 
difference in the mass flux distributions in the inviscid and 
Rey s = 10 5 cases. Also shown in the figure are the transformed shock 
layer thicknesses; 8 for the three cases. In the inviscid case; the 
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stagnation line analysis yielded a value for 5 of 0 . 796 , which is 
close to the value correlated by Inouye (o = O.jd) from inviscid 
analyses for the entire subsonic flow field surrounding a spherical 
body. In the Rey_ = 10 case, the shock layer thickness, 5 is 

o 

about 2 percent greater than that for the inviscid case, whereas for 
the Rey s = 1 case, the shock layer thickness increased about 
25 percent. 

The influence of the blowing rate on the velocity gradient and 
the shock layer thickness, y s , for a constant density is shown in 
Figure 9 where a is plotted in the physical coordinate system. The 
figure indicates that the inner flow region from the body to the stag- 
nation point (y g - (y) pv _Q & O.Okl) is drastically modified due to 
blowing but that the inviscid outer flow is virtually unaffected by 
the blowing rate. For the constant density model, the shock wave 
simply moves outward from the body as the blowing rate increases while 
the inviscid outer flow retains the same character independent of the 
blowing rate. 

The constant density results demonstrate that, for a wide range 
of blowing rates and Reynolds numbers, the solution to the flow equa- 
tions converges in a minimum number of iterations, the nodal spacing 
requirements are indeed moderate, and the computational times are 
reasonably short. With this knowledge, the behavior of the overall 
solution with the variable density iteration can now be examined. 

Non-Radiating Air Solutions 

The overall convergence of the flow solution for a variable density 
is examined in this section for an equilibrium air-gas mixture without 
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radiation. Solutions for both non -blowing and the blowing of air into 
air are presented. The air model consists of elemental mass fractions 

of O.78 nitrogen and 0.22 oxygen. 

The development of stability criteria is difficult even for the 
simplest of equations. For the coupled system of nonlinear equations 
which are solved in this study, it becomes necessary to rely on the 
behavior of the numerical solution in order to obtain infonnation of the 
overall stability and convergence behavior of the solution. The results 
of the information are reflected in the logic presented in the flow 
diagram for the overall solution procedure shown in Figure 6. It should 
be emphasized that Figure 6 represents the results of the stability 
study which is used in all subsequent studies. The numerical results 
which led to this particular iterative procedure for the variable den- 
sity solutions are discussed below. 

The solutions to be presented below are for the following free- 
stream conditions unless otherwise noted. 

Ui = 14.6 km/sec 

•p' = 1.6 x 10“^ atmos 

•^00 

p' = 2.38 x 10-7 gm/crn^ 

00 

R' = 3^2.7 cm 
b 

It was observed in the numerical solution that major oscillation 
occurred in the enthalpy and the density, particularly in the viscous 
region of the flow, as the iterative procedure progressed in time. 

These enthalpy oscillations are shown in figure 10, which is a plot of 
the calculated enthalpies at the first three nodal points adjacent to 
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the wall and within the shock layer (Note that = 0.1 and is 

constant. ) versus iteration number, where one iteration represents one 
complete pass through the governing equations. From the figure it can 
he seen that the amplitudes of the oscillations generally decrease with 
T) and, while it is not shown here, it was observed that after about 
five iterations the enthalpy values had converged beyond tj « O.J. 
However, the solutions near the wall did not converge but merely contin- 
ued to oscillate even after thirty iterations. 

The density profiles exhibited a similar oscillatory behavior as 
the enthalpy profiles but in an inverse fashion (i.e., an overprediction 
of the enthalpy led to an underprediction of the density). This is 
because there is a strong inverse coupling between the density and the 
enthalpy in the equation of state. The equation of state for air has 
been correlated by Smith, G. L. (see Garrett, Smith, and Perkins) in 
the form p « pV, where a and b are positive exponents. Since 
p is nearly constant, then an overprediction of the enthalpy obtained 
in the energy equation is a corresponding underprediction of the density 
from the equation of state and vice versa. It is possible to get into 
a resonant computing mode in which the density and the enthalpy oscil- 
late back and forth and can be slowly convergent or even divergent. 

From the results it was not clear that the solution was divergent; 
however, it was apparent that if it was converging near the wall it was 
prohibitively slow. 

Since the calculations appear to be oscillating about the solution, 
then by proper damping of the calculated quantities it is possible to 
speed convergence. Fox has discussed the method of underrelaxation 
whereby the quantity (s) which are highly oscillatory are weighted (or 
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damped) before they are used in subsequent calculations. It was found 
that by weighting the enthalpy and the density profiles by a certain 
percentage of their former value and their newly calculated values, 
rapid convergence of the overall solution could be obtained. The 
weighted values of the enthalpy and the density at each nodal point are 
computed from the following relations: 

h = (HDAMP)h ^ 1_1 ^ + (1 - HDAMP)^ 1 ^ 
n n n 

and 

p = (RODAMP)pf 1_l) + (1 - RODAMP)p^ 1 ^ 
n n n 

The convergence behavior of the solution for damping factors of 
0-5 and 0.9 is shown in Figure 11. Plotted in this figure are the 
weighted enthalpy and the density values computed at nodal point 2 . 

Both solutions converged; however, it is apparent that, while over- 
damping will insure convergence of an otherwise oscillatory solution, 
it can be unduly time consuming. For most of the cases examined, both 
with and without blowing, the 0.5 damping factors appeared to be near 
optimum in terms of generating a converged solution in a minimum number 
of iterations. For damping factors around 0.2 to 0 . 5 , certain solutions 
would converge but required more iterations than the 0.5 damping factor 
cases. It was also observed that for the largest blowing rate (-0.2) 
considered in this analysis, damping factors of 0.7 speeded convergence 
of the solution by about a factor of 2 over the 0.5 damping factor 
solution. Because of the sensitivity of the computer run times to the 
damping factors, and since no attempt was made in the present study to 
optimize the damping factors for a particular solution beyond that 
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which was just discussed., the reader is reminded that subsequent dis- 
cussions of computer time relate closely to the 0*5 damping factors. 

With a thorough study of damping factor requirements for optimum solu- 
tions, run times may be .improved significantly. 

In order to obtain an indication of the stability of the variable 
density solution for various blowing rates, three cases for blowing 
rates of 0, -0.1, and -0.2 were examined. The same free-stream condi- 
tions given previously and the 0.1 wall enthalpy value was specified. 
Twenty-one nodal points were used across the shock layer. All solutions 
converged to an absolute density accuracy of 0. 01 for damping factors of 
0.5. It was observed that the number of iterations of the governing 
equations, and consequently the run times, increased with increasing 
blowing rate. For the (pv) w = 0, -0.1, and -0.2 cases, the number of 
iterations required to obtain convergence were 11, 1 6, and 21, respec- 
tively, and the corresponding computer times were 1.5, 2.5> and. 

3.5 minutes, respectively. The computer run tunes are well within 
reason for the fully viscous shock layer with equilibrium air chemistry 
computations. 

The results obtained for the (pv) w = 0 and -0.1 cases are shown 
in Figure 12. On a physical scale the boundary layer for the (pv) w = 0 
case occupies about 5 percent of the shock layer and the combined inner 
layer and boundary layer for the (pv) w = -0.1 case occupies about 
20 percent. The (pv ) w = -0.2 results are not shown in the figure. 

The (pv) = -0.2 solution converged) however, in the inviscid outer 
region the enthalpy profiles were extremely irregular. Enthalpy values 
ranging from 0.4-5 up to about 0.7 (the latter value being greater than 
the free-stream total enthalpy) were observed. To a similar degree. 
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density profile irregularities and minor inconsistencies in the velocity 
gradient profiles were observed. Initially, it was felt that this 
oscillatory behavior of the enthalpy was due to the non-radiating gas 
assumption because the conditions were sufficient to produce non- 
adiabatic effects in the outer inviscid region. About this time in the 
study, the radiation computer code was fully incorporated into the pro- 
gram and results were being generated with the radiation flux term 
included in the energy equation. And, indeed it was observed that the 
slight oscillations which were present in the enthalpy distribution in 
the outer inviscid region for the (pv) w ~ 0 8X1(1 cases (Fig* 12 (d)) 

were not observable in the radiating flow field solutions. 

It appeared that the solutions were indicating that the radiative 
fluxes must be taken into account in the energy equation at the 
l4-55 km/sec velocity. Thus it was decided to run a lower velocity 
case of 10 km/sec where radiation is not significant, and to observe 
the behavior of non-blowing and massive blowing solutions, particularly 
the enthalpy distributions. The results of this study are shown in 
Figure 13 for (pv) = 0 and —0.2 blowing rates. The numerical 
results show that while there is improvement in the solutions for the 
lower velocity case, the oscillatory behavior of the energy equation 
solution persisted. In retrospect, the improved behavior for the lower 
velocity case is due to the decrease in the Reynolds number ( 7.8 x 10^ 
and k . 0 x 10 ^ for the = l4 . 55 W sec and the Ud, = 10. 0 Ion/ sec 
cases, respectively) which permits a coarser nodal spacing. 

It was observed much later in the study, when the solution to the 
elemental diffusion equation was required for the foreign ablation 
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product injection study, that the energy equation and the elemental 
diffusion equations are similar and have coefficients on the derivative 
terms which under further dimensional analysis reduce to about the same 
order of magnitude. As will be subsequently discussed, the elemental 
diffusion equation can be unstable unless the nodal spacing is made pro- 
hibitively small (At) « 10"^ ) and a windward difference formulation is 
required to obtain meaningful profiles. The one notable difference 

in the two equations is the appearance of the radiative heat flux and 
term in the right-hand side of the energy equation. When radiation was 
included, this term provided sufficient damping in the computations to 
yield accurate solutions to the energy equation. However, as a conse- 
quence of the elemental diffusion equation study, a windward difference 
form of the energy equation was also examined. Certain radiating cases 
were rerun, replacing the central difference form of the energy equa- 
tions, which is used in most of the results to follow, with the wind- 
ward difference form. It was noted there were no significant differ- 
ences in the radiative heat flux at the wall. Comparisons of the 
results for the two formulations appear in a subsequent section. 

Radiating Air Solutions 

The viscous radiating shock layer solutions, including air-to-air 
injection are discussed in this section. The present results are com- 
pared with results from existing approaches. 

General Results From the Present Analysis 

The results obtained for (pv) = 0, -0.1, and -0.2 cases are 
shown in Figure l4 for the following free -stream conditions: 

Uj, = 15.25 km/sec 

= 2.72 x 10~7 gm/cm^ 
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p' = 1.95 x 10“^ atm 
00 

and a nose radius of 3*048 meters. 

A wall enthalpy of 0. C28 was specified which corresponds to a wall 
temperature of 3,600° K. This temperature is close to the steady state 
ablation temperature predicted by Smith et al. (1970)* Twenty-one 
nodal points were used for all computations with the exception of the 
radiative flux computations where 11 points were used to avoid excessive 
computer time. The intermediate values of the radiative heat fluxes 
required in the energy equation solution were obtained by linear inter- 
polation. Typical computer time ranged from about 30 minutes for the 
no-blowing case to about 70 minutes for the (pv) w = -0.2 cases. 

A comparison of the radiating solutions shown in Figure l4 (a) with 
the non-radiating solutions shown in Figure 12(a) for slightly different 
free-stream and wall enthalpies shows that the physical shock displace- 
ment distance significantly decreases when radiation is taken into 
account. However, the mass flux distributions were only slightly 
altered in the T)-coordinate system. The significant decrease in y s 
when radiation is included is due largely to the non-adiabatic effects 
in the outer inviscid region which results in higher values in density 
(by alm ost a factor of 2 as shown by comparing Figures 12(e) and 14(e)) 
and to a somewhat lesser extent, the adjustment of the inner flow to the 
wall boundary conditions. On a physical scale the outer inviscid region 
occupies about 95, 75, and 60 percent of the total shock layer for 
blowing rates of 0, -0.1, and -0.2, respectively. 

The y-mamentum equation solution given in Figure l4(b) shows that 
in al i cases the pressure, as expected, reaches a maximum at the stag- 


nation point. 
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As previously mentioned, the inclusion of the radiative flux term 
in the energy equation has a damping effect on the solution. A com- 
parison of the velocity gradient profiles shown in Figures 13(a) and 

14(c) for (pv) = -0.2 shows that the waviness in the velocity pro- 
w 

file vanishes when radiation is included. Further, the enthalpy 
profiles (Fig. 14(d) for (pv) = 0 and -0.1) are regular, and although 

W 

the enthalpy profile for (pv) = -0.2 exhibits a slight waviness, it 
is considerably improved over the large oscillations shown in 
Figure 13(b) for the non-radiating cases. It is noted that this 
improvement was obtained at an even higher Reynolds number than the 
non-radiating cases (Rey_ = 9 x 10^ for the radiating case). 

The radiative heat flux distributions across the shock layer are 
shown in Figure l4(f). The positive heat flux values indicate that 
the net radiative heat transfer is away from the body and the negative 
values indicate toward the body. The nondimensional heat fluxes at 
the wall are -0.04-34, -0.0390, and -O.O369 (^ w = -^150, -37^0, and 
-3490 watts/cm^, respectively) for the (pv) w = -0.1, and -0.2 cases, 

respectively. Thus for air-to-air injection at (pv) = °> -0.1, and 0.2 
the radiative heating rates are reduced 10 and 18 percent, respectively, 
below the non -blowing rates. The absolute values of heat fluxes for the 
blowing rate cases decrease slightly in the inner region in the direction 
of the body} however, the results indicate that the inner air layer is 
only moderately effective in absorbing the incident radiation from the 
high -temperature outer layer. 

It is observed in Figure l4(f) that the radiative heat fluxes reach 
a minimum near the stagnation point, which is expected based on an 
examination of the governing differential equations. This feature 
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provides an important self -test of the overall accuracy of the present 
method, i.e. , the ability of the method to generate sufficiently 
accurate thermodynamic properties which, when input into the radiative 
heat flux computations, yield the minimum in the fluxes at the stag- 
nation point. 

More interesting and important checks on the adequacy of the 
present method for the stagnation line solution can be made by com- 
parison with solutions to the entire subsonic flow field and by com- 
parison with existing stagnation line solutions. The first exercise 
serves not only as a check on the accuracy of the solution but also as 
a check of the fundamental assumption that the stagnation line solution 
can be decoupled from the entire region of influence. 

Comparisons of Non-Blowing Results With Existing Solutions 

Shown in Figure 15 is a comparison of the enthalpy profiles 
generated in the present approach with the stagnation line enthalpy 
profiles obtained by Suttles (1969) from a one-strip method of integral 
relations solution and with those obtained by Falanga and Sullivan 
from an inverse method solution. 

The solutions axe for non-blowing, with the following free-stream 
conditions : 

IT = 14.55 km/sec 
p' = 2.377 x 10-7 gm/crn^ 
p' = 1.6 x 10"^ atm 

OO 

and a nose radius of 3-427 meters. Both of the comparison solutions 
are extracted from complete radiating solutions to the entire subsonic 
flow field. The present solution is viscous, whereas the comparison 
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cases are for inviscid flows } however, the comparisons are still meaning- 
ful since for non-blowing conditions the boundary layer occupies only a 
small percentage of the shock layer. The RATRAP computer code was used 
in all three studies. The present solution for the enthalpy profile 
compares favorably with the inverse solution in the inviscid region for 
both N = 11 and N = 21 points with radiative flux calculations at 
11 points. While the inverse solution is considered more accurate than 
the method of integral relations solution in defining shock layer pro- 
files, it is interesting to note that the radiative heat fluxes pre- 
dicted at the wall by the three methods fall within 5 percent of each 
other. 

A complete su mmar y of the radiative heat fluxes calculated at the 

wall is given in Table II for the cases examined in the present study. 

Also shown in the table are comparisons with existing sources. 

Comparisons of Air-to-Air Injection Results With Existing 
Solutions 

The viscous radiating solutions for air-to-air injection with 

(py) = -0.1 are compared in Figures 16 and 17 for the following 
w 

free-stream conditions: 

= 15.25 km/sec 

= 2.72 x 10~7 gm/cm^ 

p' = I.95 x 10“^ atm 
-^00 

and a nose radius of 5. 0^8 meters. The wall enthalpy is 0.028 which 
corresponds to a wall temperature of 5,600° K. 

Shown in Figure 16 are comparisons across the entire shock layer 
with the results of Rigdon et al. (1969) and with results supplied by 
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Wilson* in a private communication from his method presented in 1970 

* 

and hy Smith* G. L. in a private communication from the Smith et al. 
(1970) method for identical conditions as those given above. In 
general* for air-to-air injection* the flow field results from all 
approaches are in fairly good agreement although there are notable 
detailed exceptions in the flow field structure which are discussed 
below. Also* the wall radiative heat fluxes computed by Wilson* Smith 
et al . and in the present analysis are within 5 percent* whereas the 
results of Rigdon et al. are about 30 percent higher. This disagree- 
ment is attributed primarily to the differences in the radiation model 
employed in the first three approaches (RATRAP) and that employed by 
Rigdon et al. (SPECS). 

Although the radiative heat fluxes predicted by Wilson are in good 
agreement with the present method* there was a notable 10 percent 
difference in the shock layer thickness predicted by Wilson and by the 
other three approaches. The possible sources of this difference are 
associated with the x-momentum equation solution and are discussed 
below. 

Tangential Velocity Gradient . There is fairly good agreement in 
all the tangential velocity gradient results in the inviscid outer 
region for all the solutions and* in general* somewhat poorer agreement 
near the body as shown in Figure 16(a). Near the body the present 
results for the velocity gradient are close to the results obtained by 

*This author is indebted to Mr. K. H. Wilson of the Lockheed 
Aircraft Corporation and Dr. G. L. Smith of the Langley Research Center 
for generating the computer solutions for direct comparison with the 
present approach. 
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Rigdon et al. (1969)- The major differences near the body in the 
present solution and those obtained by Wilson (1970) are due to differ- 
ences in the computed inner layer and boundary layer thickness as shown 
in Figure 17(a). 

The major source of the discrepancy in Wilson's inner region 
solution and the present results is probably due to differences in the 

o o 

assumptions made for the pressure gradient coefficient, (8 p/c fee ) 

The higher the pressure coefficient, the more rapidly the oncoming 
flow sweeps around the body and, thus, the smaller the standoff 
distance. The value of the pressure gradient coefficient is not well 
defined at the present. Wilson assumes a coefficient of -3.0, whereas 
the value of -2.0 is assumed by Rigdon et al. and in the present 
analysis. Smith et al. use a value of -2.5 which has been correlated 
by Inouye on the basis of inverse flow field solutions. 

A case was run to examine this effect. The pressure gradient 
coefficient was set to -2.5 and the shock was considered to remain 
concentric. Thus, the velocity gradient behind the shock became, by 
equation (42), a s = \j p/2 = 1.12. The results indicated that the 
shock layer thickness decreased 7 percent, which agrees more closely 
with Wilson's shock layer thickness. However, most of the adjustment 
occurred in the outer inviscid region (by virtue of the higher velocity 
gradients in this region). As shown in Figure 18 there was not any 
significant change in the inner region of values of the velocity 
gradient or in the extent of the inner region. 

Two assumptions are made by Wilson which are not necessary in the 
present method. First, because of numerical stability problems, Wilson 
treats the inner region as inviscid and begins his fully viscous 
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calculations near the inner region and boundary layer interface. A case 
was run assuming that the first five nodal points, up to y/y s *** 0. (A- 
were inviscid. This was accomplished by setting p = 0 at these points. 
There was no discernible differences in the fully viscous and the 
inviscid/viscous solutions. This indicates that the inviscid inner 
region assumption is valid for the massive blowing problem. 

The second assumption made by Wilson in his solution to the govern- 
ing equations is that of a constant density viscosity product across the 
shock layer given by Pp = (Pp) • Since the value Pp within the 

"W 

boundary layer for the comparison case examined was about one-half that 
of the wall value, it was decided to examine this assumption. Two 
cases were run for (pv) = 0 -0.1 in which the density viscosity 

products appear in x-momentum and the free energy equation was set 
equal to (pp) w at each nodal point. The results of the constant Pp 
cases and the corresponding variable pp cases were nearly identical. 
This is shown in Figure 19 . The results tend to indicate that while 
viscosity is important in the boundary layer, it need not be too well 
defined in the computations since the boundary layer is actually a thin 
transition region which adjusts to the outer and inner flow regions. 

Returning to Figure 16 (a), the results from the integral solution 
of Smith for the total shock layer thickness agrees fairly well with 
the present results; however, his inviscid inner layer thickness agrees 
more closely with Wilson calculations. 

Also shown on this figure are the solution obtained by Wilson and 
Hoshizaki (1969) from their integral solution approach. The comparisons 
indicate that the method, which was sufficient for the non-blowing 
radiation studies, is not adequate for the strong blowing conditions. 



Total Enthalpy . The total enthalpy results across the shock layer 
are compared in Figure 16(b). With the exception of the differences in 
the shock displacement distances which have been noted, the present 
results compare favorably with the results of Wilson. The enthalpy 
results of Smith et al. (1970) are fairly good overall. Their method is 
based on a one-strip integral relation solution for the outer inviscid 
flow coupled with an integral solution for the boundary layer and inner 
inviscid region. It is not expected to yield the detailed flow field 
structure results obtainable from the other approaches (_i.£. , Wilson 
(1970), Rigdon et al. (1969), and the present method). 

Temperature . The temperature results obtained from the four 
approaches agree fairly well in the outer inviscid region with somewhat 
poorer agreement in the boundary layer and near the wall as shown in 
Figure 16(c). Since the radiative heat fluxes are strong functions of 
temperature, the results near the wall are compared in more detail in 
Figure 17 (b). The present results compare favorably with the results 
of Wilson. It is unfortunate that, because of differences' in the 
radiation models, no direct comparison can be made with Rigdon et al. , 
who also solves the fully coupled viscous radiating shock layer equa- 
tions. However, it is apparent from their temperature profile that the 
radiation model has a strong influence on the resulting temperature and 
the radiative heat flux predictions. 

The temperature results of Smith et al. shown here indicate that 
the temperature is nearly constant in the inviscid inner region, whereas 
all other solutions indicate an appreciable drop in temperature going 
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toward the body due to radiation exchange. This author has been 
advised by Smith that their constant temperature results were due to 
a programing error which affected the energy equation solution near 
the wall and that they now expect the temperature to increase more 
rapidly with y. 

Radiative Heat Flux . The radiative heat fluxes are shown in 
Figure 16(d). The present results indicate that the inviscid inner air 
layer is relatively ineffective in reducing the radiation to the wall 
and, as a consequence, the wall heat flux results of Wilson and of 
Smith et al. , who predict thinner inner layers, agree well with the 
present results. 

A summary of the predicted heat fluxes at the wall for the various 
air-to-air injection rates is given in Table 2-A. 

It is interesting to note that the radiative heat fluxes at the 
wall converged much more rapidly than the detailed thermodynamic proper- 
ties across the shock layer. A typical example of the sensitivity of 
the radiative heat flux at the wall on the density profile is shown in 
Figure 20. It was observed that the heat flux at the wall for the 
(pv) = -0.1 case converged within *2 percent of its final value after 

W 

only five iterations, whereas the density values in the inner region 
(t) = 0.2 and 0.4) were about 50 percent below their final value. This 
indicates that the radiative heat fluxes at the wall for air-to-air 
injection may not be as sensitive to the thermodynamic properties in 
the inner region as previously expected. If one is interested only in 
the gross quantity of the heat flux at the wall, then the 1 or 2 percent 
density convergence criteria could be relaxed considerably to provide 
substantial savings in computational time for air-to-air injection 


studies. 
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As previously noted, the air in the inner region is not too effec- 
tive in absorbing the incident radiation from the high -temperature outer 
layer. However, since certain chemical specie of ablators are strong 
absorbers and emitters, the flux to the wall may be sensitive to the 
relative amounts of these species. If this is the case, then a detailed 
definition of properties across the shock layer is desired. Thus, the 
density convergence criterion of 2 percent was retained in the subsequent 
ablation study. 

Stability Study of the Elemental Diffusion and Energy Equations 

The solution to the elemental diffusion equation is required in 
the analysis of the radiating shock layer with ablation product injec- 
tion. In extending the present method from the air-to-air injection 
analysis to the ablation products injection into air, it became evident 
that problems of numerical stability existed in the solution of the 
elemental diffusion equation. Presented in this section is the stability 
study which led to the windward difference formulation for the elemental 
diffusion equation and the subsequent reexamination of the energy equa- 
tion formulation. 

Initially, a central difference form of the elemental diffusion 
equation was used and an attempt made to solve the fully viscous 
radiating shock layer equations, including ablation product injection. 

The diffusion equation solution was unstable in the first pass through 
the equation resulting in values of Kp at certain nodal points which 
were outside the ranges of 0 <_ a j, < 1.0. The program was immediately 
aborted in the equilibrium chemistry calculations because of negative 
mass fractions of the individual species. The reason for this instabil- 
ity was investigated both analytically and numerically. 
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The properties vary too rapidly across the shock layer in the 
variable density solution to perform any meaningful analytical analysis 
of the stability of the equation. However, some insight into the 
stability of the equation and the associated nodal spacing requirements 
can be obtained from a constant density, constant diffusion coefficient 
assumption. 

For a constant P and D-]p across the shock layer, the elemental 
diffusion equation (39) becomes 

~ 2 2 ~ 
dcc.p P hq2 ^ a F 

p v = — o" 

dr] 5 di)2 


Let 



dr] 


thus. 


/ o 

p ^ d 12 

where the primes are used to denote the derivative. 

Central differencing the two above equations yields 


a n+l - a n-l _ y n+l + y n-l 
2 At] 2 


and 


7n+l - 7n-l 
2 AT] 


pvS /?n+l + 7rt-l \ 
P 2r >12\ 2 / 
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In matrix form the above equations become 


a n+l 


7n+l 



For stability, the determinant of the matrix must be positive (see, e.g. , 
C r andall ) . Thus , 


1 - 


1 + 


PvS 


p2D 


AT] 


12 


> 0 


Pv5 

p 2 d 12 


AT) 


The nodal spacing requirement for stability can be established when it 
is noted that 


6 = 0[l] and pv = 0 ['ll 

For a constant density of 20, the nodal spacing requirement becomes 


AT) < 


P 2 D i2 

ii-oo d 12 

PvS 

0 

DO 0 

0 


o[4oo d 12 ] 


Since D 12 is typically of order 10”^, then Ai)< 1) x 10“^. 



This nodal spacing requirement is prohibitive in this analysis 
since it would require at least 2,500 nodal points across the shock 
layer to generate a stable solution. However, a case was run in which 
1,001 nodal points were used for the elemental diffusion equation solu- 
tion with the same 21 points used for the other equations. Although 
some improvement was noted in the solution (for which only about 3 seconds 
was required to solve the diffusion equation), it did not eliminate the 
instability. 

The numerical results which substantiate this stability analysis 
are shown in Figure 21. The computed mass fractions for a diffusion 
coefficient of 10 - ^ is shown in Figure 2l(a) for a model spacing of 0.05* 
The solution is clearly unstable. The results obtained when the 
coefficient was increased to 1.25 x 10 and Atj was set to 0.1 (which 
accomplishes the same stability effect as decreasing the model spacing 
to 4 x 10"^) is shown in Figure 21(b). The sharp oscillations have 
vanished, but still exceeded the value of 1.0. The case was rerun 

for a nodal spacing of 0.02 with no noticeable improvement. Only when 
D -^2 was raised to 10"5 did the a-p monatonically decrease from 1. 0 
at the wall to 0.0 at the shock. 

Since the central difference form of the elemental diffusion 
equation was inadequate, a windward difference form for the convection 
term was employed. Windward differencing provides automatic damping of 
the profiles and has been shown in time -dependent studies of fluid and 
heat flow problems to be required in order to obtain meaningful solu- 
tions (see, e.£. , Richmeyer; Larkinj and Khajeh-Nouri). 

A physical explanation given for using windward differencing is 
that if fluid is flowing downstream the upstream cell influences the 
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downstream cell more than the downstream cell influences the upstream 
cell. It is correct for supersonic flow that the downstream conditions 
have no influence on the upstream flow, whereas the upstream conditions 
may have a pronounced effect on the downstream flow. Even for the stag- 
nation line where the flow is subsonic and each cell or nodal point 
influences all other cells, the explanation that the upstream cell 
influences more than the downstream cell is palatable. 

The mathematical reason for windward differencing is that higher 
order harmonics that are introduced into the solution because of finite 
difference approximations to the differential equations decay exponen- 
tially, whereas for central difference formulations these higher order 
disturbances can be exponentially amplified. 

The results obtained for both central and windward difference forms 
of the diffusion equation are shown in Figure 22. The computations were 
for a viscous radiating solution assuming binary diffusion coefficient 
based on an atomic hydrogen -atomic nitrogen mixture . The windward 
difference formulation yields a stable solution in which cc monatonic- 
ally decreases across the shock layer. 

A check case was run to establish the effect of the damping on 
the a F profiles that is introduced by the windward difference formu- 
lation. This damping, which is a form of artificial viscosity, tends 
to smooth or smear out the gradients more than is natural. To check 
this, both central and windward difference solutions were obtained for 
diffusion coefficients which yield stable solutions for central differ- 
ence approximations. The results from these two formulations, shown in 
Figure 23, indicate that the windward difference formulation does not 
introduce any noticeable artificial viscosity. 
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The enthalpy profiles in the non-radiating equilibrium air solution 
for the (pv) w = - 0.2 case were shown in Figure 13 to be unstable. 

These profiles exhibited an oscillatory behavior similar to the constant 
density profiles for the central difference approximations. As 

previously mentioned, the enthalpy profiles became smooth when the 
radiation term was included, and at that time no more attention was 
devoted to the non-radiating solutions. However, after the elemental 
diffusion equation solutions were obtained, it was decided to reexamine 
the energy equation, neglecting radiation, to see if the central differ- 
encing assumption was responsible for the oscillations in enthalpy. 

An order of magnitude analysis of the coefficients on the second 
derivative term of the energy and elemental diffusion equation revealed 
that these coefficients could be the same order of magnitude. This 
analysis is given below. 

For constant density and constant transport properties assumption, 
the energy and the elemental diffusion equations (eqs. ( 38 ) and (39)); 
respectively, can be written as (where the convective term on the right- 
hand side of the energy equation has been neglected). 


dH _ P[i d 2 H 

dq SReyPr dq2 


and 


dcc-TTi 

pv — ^ = 
dq 



,2~ 

PH ^ a F 
SReySc dq2 


where Sc, the Schmidt number, is given by 


_m7 

P'D 


12 


Sc = 
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Since the Prandtl and Schmidt numbers are typically of order unity, 
then the energy equation and the elemental diffusion equation are con- 
trolled by the same stability requirements. 

As a result of this analysis, the convective term in the energy 

^ was recast in windward difference form and the non- 
radiating equilibrium air solution for (pv) w = “0-2 was rerun. The 
enthalpy, density, and velocity gradient profiles became smooth. The 
x-momentum and energy equation solution results are plotted in Figure 2k 
and compared with the previous central difference solutions. 

A significant change in the shock layer thickness also occurred 
when the windward difference form of the energy equation was used to 
improve the profiles. The shock layer thickness, y g , was reduced by 
about 12 percent from 0.0703 to 0.0622 in going from central to windward 
difference form. It was found that this reduction in y was due 
primarily to the large blowing rate which introduced the strong 
oscillations in the profiles for the central difference formulation. 

The (pv) = 0 and (pv) = -0.1 non-radiating cases were rerun with 
w w 

the improved formulation of the energy equation. There was no change 

in y„ for (pv) = 0 and there was a 3 percent reduction in y for 
s w s 

the (pv) = -0.1 case, 
w 

The strong oscillations in the enthalpy profiles were not present 
in the radiating solutions. However, the (pv) w = -0.1 case for 
UoJ, = 15.25 km/sec was rerun to determine the effect of the windward 
difference approximation in energy equations. The radiation flux at 
the wall was unchanged and shock standoff distance increased by only 
2 percent. 


f dH 
equation Pv — 
V dT) 
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As a result of this study it was decided to retain the windward 
difference form approximation for the radiating ablation products 
analysis which follows. 

Radiating Flow Field Solutions With Ablation Products 
The solutions to the fully coupled viscous radiating flow field 
with ablation products injection are presented in this section. The 
results of the present analysis are compared with existing solutions. 

The binary diffusion model assumption is also investigated in this 
section. 

The heat-shield material is considered to be a carbon phenolic 
ablator composed of 90 to 95 percent carbon and the remaining elemental 
mass fractions consisting of nitrogen, oxygen, and hydrogen. Exact 
values of the elemental mass fractions of each of the constituents are 
specified for the particular cases presented. 

Check on the Binary Diffusion Model Assumption 

The binary diffusion coefficients were computed on the basis of 
the dominant species present. For the ablator, the dominant species 
is atomic carbon and for air it is atomic nitrogen, hence the binary 
diffusion coefficient for atomic carbon atomic nitrogen diffusion was 
used in the computations. 

However, to obtain an indication of the validity of the assumption, 
two cases were run at identical conditions with the exception that 
atomic carbon/atomic nitrogen was used in one solution in calculating 
the diffusion coefficient and atomic hydrogen/atomic nitrogen was used 
in the other solution. A comparison of the elemental diffusion equa- 
tion solution for the ablator mass fraction profiles with the two 
assumptions is shown in Figure 25 for the conditions noted there. The 
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two profiles show that the results differ in the mixing region and its 
extent, with the ablator elements extending further into the shock layer 
for the hydrogen-nitrogen diffusion coefficient. This is to be expected 
since ^ is almost an order of magnitude larger than D c _^, as shown in 
the figure. 

A detailed analysis of the heat flux distributions for the two 
cases indicated a maximum difference in the radiative heat fluxes of 
3 percent which occurred in the outer inviscid mixing region (r) = 0.7), 
as shown in Figure 26. At the wall, the nondimens ional heat fluxes 
differed by less than 1 percent (-0.02892 and -0.02917 f° r %-N 9X1(1 
D c _ n , respectively). The shock layer thicknesses, y s , were 0.04-173 
and 0.04-168 for D H _ N , and D C _ N , respectively. Thus, it appears for the 
conditions examined in this study that the binary diffusion model 
assumption .is valid. 

General Results From the Present Analysis 

The results obtained for (pv) w = 9, -0.1, and -0.2 cases are 
shown in Figure 27 for the following free-stream conditions: 

Ujjj = 15*25 km/sec 

p' = 2.72 x 10 gm/cm^ 

OO 

p' = 1.95 x 10*^ atm 

CO 

and a nose radius of 3*04-8 meters. 

A wall enthalpy of -0.04-9 was specified which corresponds to a 
wall temperature of 3,600° K for the carbon phenolic ablator injection 
cases with the following elemental mass fractions: 
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a c = 0.9207 
= 0. 0086 

a = 0.04-91 
0 

oL = 0.0216 

For the non-blowing case, air is adjacent to the wall for which an 
enthalpy of 0.028 corresponds to the wall temperature of 3,600° K. 

In general, the profiles and parameters are much the same as those 
for the air injection cases with the exceptions of the enthalpy profiles 
for the (pv) = -0.2 case and the radiative heat flux predictions. 

Even with the windward difference form of the energy equation, the 
enthalpy profile for the (pv) = -0.2 blowing rate was irregular near 
the stagnation point. This solution had not completely converged when 
the calculations were terminated on the computer after an accumulation 
of nearly 2 hours of machine time. However, it was near convergence. 
Apparently, this blowing rate is near the stability limit for obtaining 
the technique employed here for the solution to the governing equations. 
As can be seen in Figure 27(f), the density and temperature profiles, 
which are used in the radiative heat flux computations, are much 
better behaved than the enthalpy profiles. 

It can be seen from the radiative heat flux profiles of Figure 27(g) 
that the flux out of the shock (at =1.0) increased with increasing 
ablator mass injection, whereas for the air-to-air injection cases the 
flux at the shock was nearly insensitive to the blowing rate (see 
Fig. l4 -(f )). The heat fluxes at the wall, which are of primary interest, 
decreased significantly with blowing rate for the carbon phenolic 
ablator (66 and 6l percent of the non -blowing radiative heat flux rate 
for the (pv) w = -0.1 and -0.2 cases, respectively). Thus, the 
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ablation products are considerably more effective than air in reducing 
the heat flux to the wall. 

It is not clear if the large negative heat flux near the stagnation 
point for the (pv) = -0.2 case is real or if it is a result of the 
uncertainty in the thermodynamic properties calculations. However, it 
is quite likely that the effects of this point on the resulting heat 
flux at the wall is washed out in the radiation exchange between this 
point and its surrounding points. This case is compared with the solu- 
tion of Rigdon et al. ( 1969 ) in the next section. 

Comparisons of Ablation Product Injection Results With 
Existing Solutions 

Two carbon phenolic injection cases were run for comparison with 
the results of Smith et al. (1970) and Chin at a blowing rate of -O.O 76 , 
and for comparison with Rigdon et al. ( 1969 ) at a blowing rate of -0.2. 
The pertinent free-stream and wall conditions are noted in Figures 28 
and 29 for these corresponding cases. In general, the results for the 
thermodynamic and flow properties are in reasonable agreement, but 
there are noticeable exceptions relating to the heat flux computations 
which are discussed below. 

In the Chin comparison case, it should be noted that the close 
agreement of the wall heat flux prediction with Chin's computations may 
be fortuitous. Chin's earlier radiation model is sufficiently different 
from RATRAP to make any quantitative comparison meaningless. This 
is evident from Table 2 where it is observed that the non-blowing 
radiation heat flux prediction of Chin and of the present method differ 
by 20 percent. 

The comparison of interest for this case is with the results of 
Smith et al. ( 197 O). It is noted in Figure 28(d) that the present 
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solution predicts a wall heat flux about 30 percent lower than that of 
Smith. An important conclusion made in their paper is that the radia- 
tion to the wall is not attenuated as much as previously predicted by 
Chin and by Rigdon et al. (1969)- This was attributed to the presence 
of large percentages of CN in the mixing region which is a strong emitter. 
A comparison of the mole fractions of the major chemical species is shown 
in Figure 28(d). The comparisons are good considering the differences 
in the temperature distributions; however, it is noted that the peak 
mole fraction for CN predicted by the present method is only about one- 
half that of the comparison case. This discrepancy in the level of CN 
occurs in the region where Smith et al. ( 1970 ) employ a cubic fit to the 
elemental diffusion equation in order to join the wall products in the 
inner shear layer with the air composition beyond the boundary layer. 

It thus appears that contrary to the air-to-air injection study results 
which did not require too accurate a resolution of properties near the 
wall, the ablation products study indicates that these properties must 
be well defined near the wall to generate accurate wall radiative heat 
flux predictions. 

The comparison with Rigdon et al. ( 19 69) indicates reasonably good 
agreement in the radiative heat flux profiles, although there are sig- 
nificant differences in the radiation models employed. The predicted 
shock layer thicknesses are within about 5 percent. A slight waviness 
is noted in Rigdon' s Pv profile (Fig. 29(a)) near the stagnation 
point (y/y s =0.3)* In all the cases calculated by the present method, 
the Pv profile was smooth in this region. It is not apparent if this 
irregularity predicted by Rigdon et al. is real or if it related to 
their numerical procedure for integrating out in both directions from 
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the stagnation point. It is noted that their solution also exhibits an 
irregularity in the temperature profile near the stagnation point as does 
the present solution. 

Summary of the Wall Radiative Heat Flux Predictions 

A summary of the wall radiative heat fluxes calculated by the 
present method for the various conditions examined in this study and 
corresponding comparison results from previous investigations are tabu- 
lated in Table 2 . With the exception of the non-blowing results of 
Suttles (1969) and the air injection comparison case with Wilson (1970), 
the present method predicts slightly to much lower heating rates than 
calculated by previous investigators. 

These results are also summarized in Figure 30, which illustrates 
the effectiveness of the carbon phenolic ablator products in reducing 
the wall radiative heat fluxes. As can be seen from the figure, there 
is general agreement of all sources that air injection is moderately 
effective in reducing the heat transfer to the wall and, with the 
exception of the Smith et al. (1970) results, that the ablation products 
of the carbon phenolic heat shield is highly effective in reducing the 
fluxes. 

As can also be seen from Figure 30 and Table 2 , there is no con- 
sistent agreement from any sources when both the non-blowing and blowing 
results are considered. For sources which use different radiation 
transport models (Rigdon et al. ( 1969 ) and Chin), this disagreement is 
due, at least in part, to the differences in the radiation model. For 
sources which use the same model. Smith et al. (1970) and the present 
method, this disagreement is attributed to the numerical procedures 
employed. There was insufficient data available to fully evaluate 
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the Wilson (1970 ) results, with the exception of the assessment of the 
validity of the assumptions required in his analysis. These assump- 
tions have been previously discussed and shown to be valid for the 
(pv ) = - 0.1 air-to-air injection case. 

For air-to-air injection, the present results indicate a 10 and 
17 percent reduction in the wall radiative heat fluxes for (pv) w - - 0.1 
and - 0 . 2 , respectively. For carbon phenolic ablation products injection 
at corresponding blowing rates, the reduction is 3 6 and 39 percent, 
respectively. Since ablation rates are typically expected to be about 
0.1 (see Smith et al. ( 1970 )), then the 36 percent reduction in the 
radiative heat fluxes represents a significant savings in heat-shield 
weight. 
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SUMMARY AND CONCLUSIONS 


An implicit finite difference scheme is developed for the fully 
coupled solution of the viscous radiating stagnation line equations 
including strong blowing. Solutions are presented for both air-to-air 
injection and ablation products injection with blowing rates up to 
20 percent of free-stream mass flow rates. The free-stream conditions 
examined are typical of interplanetary return conditions into earth's 
atmosphere near the point of peak radiative heating in the entry trajec- 
tory. A detailed radiative transport computer code (RATRAP) which 
accounts for both continuum and line radiation exchange processes is 
utilized in the study. 

Starting with a minimum number of assumptions for the initially 
unknown parameters and profile distributions, convergent solutions to 
the full stagnation line equations are rapidly obtained by a method 
of successive approximations. No singularities exist in this formu- 
lation of the finite difference equations. Damping of selected pro- 
files is required to aid convergence of the massive blowing casesj 
however, even for these cases, no patching of the viscous and inviscid 
regions is required. The results demonstrate that windward differencing 
of the convective term in the elemental diffusion equation is required 
for a stable solution to this equation. While the central difference 
approximation to the energy equation yields satisfactory solutions 
when radiation is included, the results are considerably improved for 
blowing rate of 0.2 when the convective term is windward differenced 
in this equation also. 

Comparisons are made with currently existing solutions to the 
radiating shock layer problem. The present method predicts lower 
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wall radiative heat fluxes for carbon phenolic ablation than those pre- 
dicted by previous investigators. 

The results indicate that the ablation products are highly effec- 
tive in blocking the incident radiation from the high-temperature outer 
layer of the shock. For blowing rates of 0.1 and 0.2, typical reductions 
range from about 35 to 1+0 percent of the non-blowing radiative heat 
fluxes at the wall. 

The inner air layer is shown to be relatively ineffective in block- 
ing the incident radiation, hence the thermodynamic properties need not 
be as well defined for air injection as for ablation product injection. 

The binary diffusion model assumption was examined in the present 
analysis and, while a multicomponent diffusion study may remain meaning- 
ful for conditions and ablators other than those examined in the present 
study, the results indicate that the multicomponent diffusion model is 
not required. 

The present results are sufficiently encouraging to recommend that 
the present method be extended to radiation calculations in the presence 
of chemical nonequilibrium. Since the solution to the governing flow 
equations requires only about 5 percent of the total computational 
time, with the radiation flux computation comprising the remaining 
95 percent of the time, the implicit finite difference scheme should be 
relatively efficient for performing the computations. 


\ 

Table 1. Diffusion coefficient 
molecular constants 


Molecular constants 

Species 

e/k, °K 

0 , X 

M, gm/gm mole 

C 

50.6 

3.385 

12.01 

N 

71 A 

3.298 

A. 01 

H 

37-0 

2.708 

1.008 
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Table 2. Summary of radiative heat fluxes at the wall 
A. Air Injection 

1. Ui = 14.55 km/sec, Pi = 2-577 x 10" 7 gm/cm^, pi = 1.6 X 10~^ atm, 
R^ = 5-427 m 


(PV)„ 

qJ (watts/cm 2 ) 

K j W 

Suttles (1969) 

Falanga and Sullivan 

Present Analysis 

0 

-2860 

-5000 

-2950 


2. Ui = 15-25 km/sec, Pi = 2-72 x 10' 7 gm/cm3, pi = 1-95 X 10“^ atm, 
R^ = 5-048 m 


( pv )w 

q' „ (watts/cm 2 ) 

Kj W 

Rigdon et al. (1969 ) 

Smith et al. (1970) 

Wilson (1970) 

Present 

analysis 

0 

-5990 

-4170 



-4150 

-0.1 

-5150 

-5991* 

-5660* 

-5740 

-0.2 

-4650 


— — “ 

-5490 


B. Carbon Phenolic Injection 

1. Ui = 15.25 km/sec, Pi = 1-77 X 10" 7 gm/cm3, pi = 1.21 X 10"^ atm, 
R^ = 2.56 m 


MW 


q R,W (w^W 0 ™ 2 ) 

Chin 

Smith et al. (1970) 

Present analysis 

0 

-5050 

-2795 

-2540 

-0. 076 

-1642 

-2188 

-1620 


2. U' = 15.25 km/sec, Pi = 2.72 x 10'T gm/cm3, pi = I .95 x 10"^ atm, 
R^ = 5.048 m 



q R,W (watts/cm 2 ) 

{P V 'w 





Rigdon et al. (1969) 

Smith et al. (1970) 

Present analysis 

0 

-5990 

-4170 

-4150 

-0.1 

— 

— 

-2790 

-0.2 

-2670 

-5246 

-2560 


♦Private communications. 
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Figure 1 . Diagram of various assumptions for the analysis of radiating 
shock layers (from Anderson, 1968, p. 2 ) 
















Node number N 


Notes 



: (l) N must be odd 

(2) Nodal points are equally spaced in q 


Figure 5 - Finite difference network 




Figure 6. Flow diagram of the overall solution procedure 
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Figure 8 


(a) Velocity gradient distributions 

Constant density solutions for inviscid and viscous flows 
without blowing 
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(b) Mass flux distributions 
Figure 8. Concluded 
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Figure 11. Effect of profile dam; 

variable density solu - 
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(b) Density 
Figure 11. Concluded 
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(b) Energy equation solution 


Figure 13. Concluded 




(a) Continuity equation 

Figure l4. Equilibrium air solutions with radiation included 
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(b) Y -moment-urn equation 
Figure l4. Continued 
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(d) Energy equation 
Figure l4. Continued 
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(e) Equation of state 
Figure l4. Continued 
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(f) Radiative heat flux 
Figure 14. Concluded 
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Figure 15. Comparison of the enthalpy profile along the stagnation line with 
solutions from complete subsonic flow field calculations - 
equilibrium air with radiation and without blowing 
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(a) Tangential velocity gradient 

Figure l6. Comparison of equilibrium air solutions including radiation with 
air to air injection 
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Figure l6. 

Continued 
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(a) Tangential velocity gradient 

Figure 17- Detailed comparisons near the wall of existing equilibrium air 
radiations shock layer solutions with air to air injection 
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Temperature, 



( b ) Temperature 


Figure 17- Concluded 
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Figure l8. Effect of the pressure gradient on the tangential velocity- 
gradient solution 
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(a) Tangential velocity gradient 


Figure 19* Effect of constant density viscosity product assumption on the solution 
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(a) Density behavior 

Figure 20. Study of the sensitivity of the . radiative heat flux to the 
density as the solution converges - equilibrium air with 
radiation 
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(b) Radiative heat flux behavior 


Figure 20. Concluded 
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22. Stability of the elemental diffusion eqi 
and windward difference approximations 




Figure 23. Comparison of solutions of the elemental diffusion 
equation for central and windward difference 
approximations in the region where the equations 
are stable 






(a) X-momentum equation 


Figure 2k. 


Effect of windward differencing the energy equation on the 
non-radiating air solution with massive blowing, 

'pv) w = -0.2 
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(a) Continuity equat 


jure 27- Radiating shock layer solution: 
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(t) Y -momentum equation 


Figure 27- Continued 
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(d) Energy equation 


Figure 27- Continued 
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(e) Elemental diffusion equation 


Figure 27- Continued 
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(g) Radiative heat flux 


Figure 27- Concluded 




15.25 km/sec 
1.77X10" 7 gm/cm 3 
1.21X10' 4 atm 
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O Present method 
□ Rigdon, et. al. (1969) 
A Smith, et. al. (1970) 


V) w = o 


Q Wilson (1970) 

Open symbols are for air 
injection 

Closed symbols are for 

carbon phenolic injection 

Flagged symbols are heat 
fluxes normalized with 
the non-blowing heat fluxes 
calculated from the present 
method 

Note: Pertinent free stream conditions 
and nose radius values for 
each data point are given in 
Table n. 


Figure 30. Summary of wall radiative heat flux predictions for air 
and for ablation product injection 
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APPENDIX 
RADIATION MODEL 

The RATRAP radiative transport model developed by Wilson ( 1967 ) 
has been discussed in great detail by Wilson (1967b Suttles (1968), 
and Wilson and Hoshizaki (1969)- A summary of Suttle’s discussion on 
the development of the governing radiative transport equation employed 
by RATRAP and some details of the radiation model are given in this 
appendix. 

The RATRAP code, which includes the detailed line radiative calcu- 
lations (sometimes referred to as RATRAP II, since it has been upgraded 
by its originators, see Wilson and Hoshizaki (1969)), I s written for 
the calculation of the radiative heat flux at any point within a planar 
(tangent) slab in which the thermodynamic properties vary in only the 
direction normal to the slab. Local thermodynamic and chemical equi- 
librium is assumed. A distribution of two thermodynamic variables 
plus the elemental mass fractions for carbon, nitrogen, oxygen, and 
hydrogen gaseous mixtures is required for the radiation computations. 

For the tangent slab approximation, there is no component of heat 
flux in the x-direction. Thus the magnitude of the radiative heat flux 
vector is the radiative heat flux in the y-direction which is expressed 
as 

p 00 p 

M " q H,y =7 0 Jo M5s 5y)dfl ^ (A ' 1> 

where is the specific intensity, e g is the unit vector in the 

arbitrary direction in which I v is evaluated as shown in Sketch A, 
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and and v are the solid angle and the frequency of radiation, 
respectively. (Note that the quantities are dimensional; however, the 
primed symbol used to denote dimensional quantities in the main body of 
this study is dropped in this appendix for simplicity. ) 



surface 


Sketch A 

The differential equation governing the radiative transfer (see, 
e.g. , Vincenti and Kruger) which is given by 


dI v , 

— = a v (B v - I v ) 


(A-2) 


has the solution 


I v (s) = / a, v (|)B v (|)e ^ d£ + i v (o)e 

d n 


v ( e )d e 


a. 


,U)dt 


(A-5) 


where a v is the modified linear absorption coefficient, B y is the 
Planck function given by 



l6l 


_7 -Kv^ (-ftv/kT) 1 1 ), >, 

B s 2 X 10 ' -3- e - 1 (A-4 ) 

v L 

I (0) is the radiation intensity coming from the wall and e and I 
are dummy variables for the integrals such that £ < 6 < s and 
0 < | < s } as shown in Sketch A. 

Equations (A-l) and (A-3) can be combined (see Suttles (1968)) to 

yield 



when the heat fluxes into the slab at the boundaries are neglected. 


The quantity Eg 


f ! de ) 


is the exponential integral described by 


Kourganoff and is given by 



In order to evaluate equation (A-5)> it is necessary to calculate 
the spectral linear absorption coefficients a v which are functions of 
the thermodynami c properties (p and T) and the number densities of the 
chemical species. Twenty chemical species are considered in the 
thermodynamic calculations. They are C 2 > 3^2^ 0, H, CO, 

CN, C 2 H, CjH, C^H, HCN, C 2 H 2 , C - , C + , N + , 0 + , and H + . 

The total spectral absorption coefficient is separated into 
continuum and line contributions in the RATRAP calculations. The 
continuum spectral absorption processes considered are the free-free 



transitions (acceleration of free electrons in the vicinity of 
atoms and ions) of C, C , N, N , 0, 0 and the free-bound and 
hound-free transitions (deionization and ionization, respectively, 
of the H atoms and positively ionized particles. Molecular 
bands are also included in the continuum calculations. These 
molecule contributions and frequency intervals from Wilson and 
Hoshizaki (1969) are tabulated below. 


Table A. Molecular band systems used in RATRAP 


Molecule 

contribution 

Frequency range for significant 
absorption photon energy, eV 

Hg Werner 

H < «v < 15.494 

Photoionization 

15.494 < ftv < 25 

Cg Swan 

1.8 < 'ft v <6.0 

Fox Herzberg 

1.8 < ftv < 5.35 

Mulliken 

5.35 < ftv < 6.0 

Freymark 

i.8 < ftv < 6.0 

CN Violet 

2.0 < ftv < 6.0 

CO 4 th Positive 

7 < ftv < 10 

N 2 Birge-Hopfield 

11 < ftv < 14.2 

0 o Schuman-Runge 

7 < ftv < 9.2 




For the calculation of the continuum contribution to the heat 
flux, the absorption coefficients for the individual species are 
weighted with their respective number densities and approximated by 
curve fits over frequency ranges. The heat flux equation (A- 5 ) is 
subsequently evaluated by numerically integrating over y and over 
frequency, using 11 values of y and 31 values of frequency. 
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The atomic line radiation component arises as a result of the many 
bound electronic transitions which occur in the atomic nitrogen and 
oxygen species (line radiation from carbon and hydrogen species are not 
included in RA.TRAP )'. This component is obtained by first grouping 
certain line contributions within various frequency intervals. The 
net line radiation is then calculated by summing the contributions 
from these line groups. Eighteen line groups, with a total of 65 lines, 
are used in RATRAP. 

The total radiative heat flux at each point is obtained by adding 
the continuum and the line radiation contributions. 



